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ABSTRACT 

A boundary integral formulation for the analysis of 
circular plate bending under lateral loads is developed using 
Green’s functions. The formulation specifically applies to 
annular plates with arbitrary boundary conditions. The plate 
bending solution in the plastic range is determined using a 
numerical method of incremental loading. A computer program to 
perform the required calculations was developed and is 
presented. Results for three case studies are included and 
compared with results obtained by other methods. Plate 
behavior in the elastic range is in excellent agreement with 
other analytical solutions, and in the plastic range is in 
reasonable agreement with published results obtained using a 
finite element method. 
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NOMENCLATURE 



a - 


plate outer radius 


b 


plate inner radius 


c 


plate outer radius for Green’s function 


D 


flexural rigidity 


d 


plate inner radius for Green’s function 


E 


Young’s modulus 


h 


plate thickness 


I 


area moment of inertia of beam cross section 


Kr 


radial curvature 


Kq 


tangential curvature 


L 


beam length 


M,Mg 


elastic moments for beam bending problem 


Mr>MGr 


elastic radial moments per unit length 


MpP 


plastic radial moment per unit length 


Mp 


actual radial moment per unit length 


M0>Mg0 - 


elastic tangential moments per unit length 


M0P 


plastic tangential moment per unit length 


Q.Qg 


shears for beam bending problem 


Q r > ®Gr ~ 


elastic radial shears per unit length 


Qr 


actual radial shear per unit length 


q, qG 


distributed loads per unit area 
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Rr 

Re 



Oy 

w, WG 

X 

2 

6 £ T\ 



£rP 



«8 






£gi 



£y 

i,iG 

T 

6 

K, Kq 
PG 

ffr 

<^Q 

V 

a 



radial radius of curvature 
tangential radius of curvature 
radial distance from plate center 
radial position of ring load 
yield stress 
deflections 
distance along beam 

lateral distance from plate center 

plastic strain increment 

total radial strain 

elastic radial strain 

plastic radial strain 

total tangential strain 

elastic tangential strain 

plastic tangential strain 

yield strain 

slopes 

plate boundary 

angular position in tangential direction 
biharmonic operator 

loads per unit length for beam bending problem 

magnitude of Green’s function ring load 

equivalent stress 

radial stress 

tangential stress 

Poisson’s ratio 

plate domain 
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CHAPTER 1 



INTRODUCTION 



The solution of plate bending problems in the plastic 
range has been approached using finite element, and more 
recently, boundary integral methods (also called boundary 
element methods). Unlike the finite element method, which 
discretizes the inside of the plate into a number of small 
elements, the boundary integral method discretizes only the 
plate boundary, thereby reducing the order of the problem by 
one . 



Many different approaches have been used to formulate 
plate bending problems using boundary integrals, with perhaps 
the most attractive proposed by Stern^^^ and Bezine*'2>, These 
authors, working independently, developed a direct boundary 
integal method using Green’s functions to solve the general 
plate bending problem in the elastic range. This work has been 
extended by Moshaiov and Vorus*'^^ to include plastic behavior. 

Symmetry is commonly used to reduce the size of the system 
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of equations needed to analyze a symmetric structure. In the 
case of axisymmetr ically loaded circular plates the problem 
becomes one-dimensional. For the linear elastic case, it has a 
closed form analytic solution. However, some authors have 
instead treated this circular plate problem as a two- 
dimensional problem, perhaps for the sake of demonstration 
(Kamiya and Sawki^'^^). They have used symmetry to allow 
discretization of a portion of the boundary instead of the 
entire boundary. 

This thesis develops a general formulation for solving 
annular plate bending problems using boundary integrals that 
reduces the problem to essentially a one-dimensional problem by 
the use of a "ring” type Green's function. Such a formulation 
is most useful, inspite of the closed form solution, as it 
allows the analysis of plates with different boundary 
conditions and loads with one algorithm. Moreover, it permits 
a solution in the plastic range, which is not possible with a 
conventional analytic approach. In addition, the simplicity of 
the one-dimensional formulation given here has a unique 
significance for the teaching of boundary element methods. 

Butterf ield^^^ has developed a boundary element 
formulation for the one-dimensional elastic beam bending 
problem. Here, a similar approach is taken for the annular 
plate. It is also extended to include non-linear behavior 
using an incremental load method as outlined by Moshaiov and 
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Vorus. Numerical results for three different annular plate 
configurations are presented and compared to results obtained 



using other methods. A discussion of these results follows, 
well as suggestions for future work in this area. 



as 
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CHAPTER 2 



DEVELOPMENT OF THE GOVERNING DIFFERENTIAL EQUATION 



This chapter reviews the derivation of the governing 
differential equation for a circular plate subjected to a 
symmetrically distributed lateral load. The derivation is 
based on that provided by Timoshenko^ , using Kirchoff's 
assumptions. It is, therefore, applicable only to the linear- 
elastic case. Treatment of non-linear behavior is addressed in 
Chapter 5. 



2 . 1 Moment-Curvature Relationships 

The first step in developing the governing differential 
equation for a circular plate is to find a relationship between 
moment and curvature. Curvature in the radial direction (Kj.) 
and the tangential direction (Kg) for a symmetrically loaded 
circular plate is found using geometrical considerations. , 

which is the inverse of the radius of curvature in the radial 
direction (Rj.)> can vary as a function of the distance (r) from 
the plate center. Examining the plate of Figure 2-1 at a small 
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radial distance dr away from r, the slope (♦) can be expressed 
in terms of the radial distance r and the deflection (w) as 
follows 



i 



dw 

dr 



The radial radius of curvature is 



R 



r 



dr 

d^ 



( 2 - 1 - 1 ) 



giving an expression for the radial curvature in terms of w 



1 d» d’« 

K = - = — j (2-1-2) 

R dr dr 

r 
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Figure 2-1. Deflection, slope, and curvature relationships 
in circular plate bending 



Due to symmetry, the tangential radius of curvature is 
constant at any given radius r, and for small deflections can 
be expressed as 



R 



e 




giving an expression for the tangential curvature 
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1 ♦ 1 dw 

^ Rg r r dr 

To relate moment to curvature, consider the small plate 
element illustrated in Figure 2-2. Define Mj. and Mg, the 
radial and tangential moments per unit length respectively, as 
follows : 



M 



r 



,h/2 



a z dz 
r 

■'-h/2 






h/2 

ff z dz 
-h/2 



(2-1-3) 



where h is the plate thickness. Assuming the plate is thin and 
hence experiences a two-dimensional stress state, Hooke’s law 
yields the following expressions relating stresses to strains 



<7 

r 



E 

d-v') 




<7 



6 



E 

d-v') 




(2-1-4) 



where E is Young’s modulus, v is Poisson’s ratio, and Cj.® and 
£Q® are the elastic radial and tangential strains, respec- 
tively. In Chapter 5, which discusses non-linear behavior, a 
distinction will be made between the elastic strain and the 
total strain. In the plastic range, the total strain includes 
elastic and plastic components, in accordance with the 
following expressions: 
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where the superscript p refers to the plastic strain 
components. Since Chapter 2 treats only elastic behavior, the 
total strain is simply equal to the elastic strain. 




Substituting Eq. (2-1-4) into Eq. (2-1-3) gives 
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E 



h/2 



M 

r 



M 



8 



(1-v ) 



z dz 



-h/2 



E 

d-v') 



h/2 

[ Cq + V ] 2 dz 

-h/2 



(2-1-5) 



Finally, from geometrical considerations, strain and curvature 
are related by 



e 

r 





£ 



e 






1 dw 
r dr 



( 2 - 1 - 6 ) 



which, when substituted into Eq. (2-1-5), give the desired 
moment-curvature relationship for a circular plate 



M 



- D 



. 2 

d w 



dr 



V dw 
r dr 



M. 



- D 



1 dw 



r dr 



+ V 



d ^w 



dr 



where D, the flexural rigidity, is given by 



E h 



12 (1-v ) 



(2-1-7) 



D 



( 2 - 1 - 8 ) 



2 . 2 Equilibrium Equations 

For the moment-curvature relationships to be useful, an 
equilibrium equations in terms of moments and shears must be 
found. The governing differential equation is then developed 
by substituting the expressions relating moment to curvature 
into the equilibrium equation. 

To find an equilibrium equation, consider first the 
laterally loaded circular plate element illustrated in Figure 
2-3, where Qj- is the radial shearing force per unit length and 
q is the distributed load (force per unit area). The couples 
acting along edges ab and cd due to radial bending moments are 



ab : 



M j. r d0 



dM 





cd: 



M + 



r 



r + dr 



r 



dr 



The shear forces acting along the edges are 



ab : 



Q r d8 



dQ 




cd : 




r + 



dr 



d8 




Figure 2-3. Shear forces and moments acting on differentia^ 
plate element. 



Because of symmetryj there are no shear forces acting 
along edges ad and bC| leaving only the tangential bending 
moments 



These couples can be resolved into components acting 
perpendicular to the axis ro (which are equal and opposite and 
thus cancel), and components acting along the ro axis which 
have a total magnitude of 
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M dr d9 

D 



Summing up all the couples and neglecting the small change in 
shear force across the element (i.e. dropping the dQ terms) 
gives the following equilibrium equation 



dM 

Mr- + — — dr 

dr 



r + dr 



de 



M r de 
r 



— Mg dr dS + Q r d9 dr = 0 

This equation is further simplified by ignoring higher order 
terms to give 

dM 

M + r -M. + Qr = 0 (2-2-1) 

r , o IT 

dr 



Also, from equilibrium in the z direction 



dQ 



dr 



( 2 - 2 - 2 ) 



2 . 3 Differential Equation for Circular Plate Bending 

Having found the desired equilibrium equations, the moment 
curvature relationships of Section 2.1 are now used to develop 
the governing differential equation for the circular plate 
bending problem. Specifically, substituting Eq. (2-1-7) into 
Eq. (2-2-1) results in a third order differential equation to 
describe the response of a circular plate to lateral loading 
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1 dw 



(2-3-1) 



. 3 

d w 



dr 



1 d w 
+ — 



Q 



r dr 



r dr 



which can be rewritten to provide an expression for shear 



Q 



, 3 

d w 



dr 



1 d w 

+ — 

r dr ^ 



1 dw 
r^ dr 



(2-3-2) 



Differentiating Eq. (2-3-1) with respect to r and using Eq. (2- 
2-2) results in a fourth order differential equation in terms 
of r, w, D, and q 



j 4 

d w 



dr 



2 d w 
r dr^ 



1 d w 



1 dw 



r dr 



+ — 



r dr 



q 

D 



(2-3-3) 



which can be expressed in a more concise form using the 
biharmonic operator as 



7*w 



q 

D 



(2-3-4) 



where 



dr 



2 d 



+ — 



r dr 



1 d 



2 , 2 
r ar 



1 d 



+ — 



(2-3-5) 



r dr 
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CHAPTER 3 



A BOUNDARY INTEGRAL METHOD 



This chapter describes a boundary integral method, and 
how it can be used to provide a general method of solution of 
circular plate bending problems. A one-dimensional beam 
bending problem is used to illustrate this method. 

3 . 1 General Formulation 

Solution of Eq. (2-3-3) analytically as a boundary value 
problem by applying known boundary conditions is possible 
(Timoshenko). However, this approach requires special 
attention for each plate configuration and load distribution. 
This difficulty is compounded when considering problems 
involving plastic behavior which should be solved numerically 
due to the non-linear behavior in the plastic range. 

Therefore, a general method is sought wherein the same 
procedure can be used to solve a variety of problems, thus 
lending the problems to computer solution techniques. One such 
method, adopted for this work to examine circular plate 
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bending, is the boundary integral method. 



Stern and Bezine, working independently, obtained a 
general formulation for plate bending problems in terms of 
boundary integrals. These integrals involve displacements, 
slopes, bending moments, and shears on the plate boundary. 

Both authors used the following reciprocal identity, which can 
be derived using Green’s identity: 



where 



( ”c 


4 

7 w 

t 


- 


w V^Wp 1 dfl = 


Q 


( 


W 


* Mg ♦ - M - Wg Q ) dr 




r 








Q 


- 


plate domain 




r 


- 


plate boundary 




D 


- 


flexural rigidity 






- 


biharmonic operator 




W, WQ 


- 


deflections 






- 


slopes 




M,Mq 


- 


bending moments per unit length 




Q.Qg 


- 


shear forces per unit length 



To arrive at the boundary integral equation, wq is selected as 
a Green’s function. 



3 . 2 Beam Bending Problem 



Butterfield has presented a reduced form of Eq. (3-1-1) 
for the one-dimensional case of beam bending, which is 
developed in the discussion that follows. As a first step, 
Butterfield multiplies both sides of the governing differential 
equation for a beam by a second function and integrates over 
the beam’s length, giving 



L 

El 7^w dx 

VJ 

0 



L 

\ w_ dx 

Cj 

0 



(3-2-1) 



where \ is the distributed load per unit length, I is the area 
moment of inertia of the beam cross section, and the operator 
7“* is defined for the beam case as 



4 

7 




Integrating the left hand side 
several times, Butterfield develops 



of Eq. (3-2-1) by parts 
the equation 



L 





' \ Wq 


- 


w El 7^Wg 


1 dx 


= 


0 













-Wq Q + M - ^ 




1 ^ 
•*0 



Recognizing that for the second function 

E I 7^^ = 

the preceding equation can be rewritten as 
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L 



( ^ ” 



w \ 



G 




0 



L 



-“a ® « 



Mg » + 9g w 



(3-2-2) 



0 



To gain some physical insight into Eq. (3-2-2) it is 
useful to recall Betties reciprocal work theorem. This theorem 
states that for an elastic structure subjected to two 
independent causes (e.g. loads), the total work done by the 
first cause in moving through the displacements resulting from 
the second cause is equal to the total work done by the second 
cause in moving through the displacements produced by the first 
cause^"^^. The theorem, as it applies to our beam bending 
problem, can be better understood by examining the examples of 
Figure 3-1. 
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Figure 3-1. Examples of Betti’s reciprocal work theorem 

In Figure 3-la, the product of load and displacement 
is equal to the product of load Pj and displacement Sji- 
Similarly, in Figure 3-lb, the product of moment and slope 
♦ij is equal to the product of moment Mj and slope ♦ji- 

This concept can be applied to understand Eq. (3-2-2), 
with the i cause the actual loading or applied moment condition 
for the system of interest and the j cause the loading or 
applied moment condition of the second function. In this 
equation, the K term represents the external load in the system 
of interest, which when multiplied by the deflection wq 
described by the second function, gives a reciprocal work term. 
Similarly, the 4(j term is a slope for the second function, 
which when multiplied by the moment M for the system of 
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interest yields another work term. This concept applies to all 
the terms of the equation, which are combined according to 
Betti’s theorm of reciprocal work to yield the given equality. 

For Eq. (3-2-2) to be useful for solving the beam bending 
problem, Butterfield chooses as the second function a Green’s 
function which is based on a singular loading condition; that 
is to say, a point load. If the point load of unit magnitude 
is applied at either x=0 or x=L, the second intergral term on 
the left hand side of Eq. (3-2-2) takes on the value of the 
boundary condition of the beam of interest. Usually, some 
attention should be given to cases where the integrals become 
singular . 

This yields two equations. However, in general, beam 
bending problems involve four unknown boundary conditions. For 
example, if simple supports are specified on both ends, the 
moments and deflections at the two ends are known to be zero, 
while the shears and slopes are not known. Therefore, two 
additional equations are required. 

To obtain two more equations, derivatives of the first two 
equations are taken. This gives an additional Green’s 
function. Boundary conditions can now be obtained by solving 
the four equations simultaneously. Knowing the boundary 
conditions, the beam problem can be solved for interior points 
using Eq. (3-2-2) and appropriate derivatives with the point 
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load positioned at the location of interest. 



The formulation for the beam bending problem outlined 
above is a direct approach based on Green’s identity. 
Butterfield also develops a boundary integral equation using 
the more intuitive indirect method. For the circular plate 
formulation that follows in Chapter 4, only the direct method 



is used. 



CHAPTER 4 



DEVELOPMENT OF THE FORMULATION FOR CIRCULAR PLATES 





This chapter 
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Figure 4-1 depicts a symmet 


r i ca 


lly load 


ed 


annular plate 


an 


d includes notat 


ion that will 


be r 


ef erred 


to 


throughout th 



development of the integral equation and the discussion of the 
selected Green’s functions. 



Important symbols are defined as follows: 



W, Wq 

Mr > Mgr 
Ms, M q0 
Q r » ®Gr 
Qs, QG8 

a 

b 

h 

r 

q 



deflections 

slopes 

radial bending moments per unit length 

tangential bending moments per unit length 

radial shear forces per unit length 

tangential shear forces per unit length 

plate outer radius 

plate inner radius 

plate thickness 

radius of interest 

distributed load per unit area 



Figure 
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4 . 2 Development of Integral Equation 



Though the laterally loaded circular plate bending problem 
is two-dimensional, symmetry reduces the problem to essentially 
a one-dimensional type problem. Therefore, it is expected that 
the integral equation developed by Butterfield for the beam 
case can be adapted to the case of the circular plate. This 
equation (Eq. (3-2-2)) is repeated below: 






1 


[ \ w_ 




w \_ ] 


1 

n 


i G 




G ) 



0 



-"g ® * »g ” 




L 



+ ®G ” 



'0 



To obtain some indication of what the circular plate 
integral equation will look like, a non-rigorous method is 
first used to adapt Eq. (3-2-2) to the circular plate case, 
followed by a more rigorous mathematical development. 



In the non-rigorous approach, a first step is to replace 
the integral on the left hand side with an area integral and 
the line load per unit length \ with the load per unit area q. 
Next, adjustments are made to reciprocal work terras on the 
boundaries (right hand side of Eq. (3-2-2)). Specifically, the 
shear terms are replaced with the product of the radial shear 
per unit length and the total length (2sr), observing that for 
symmetrical loading, the tangential shear Q 0 is 0. Moment 
terms are replaced with radial moments, since there are no 
tangential moments on the boundaries. These moments must also 
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be multiplied by the length (2icr) to give total moment. The 
resulting expression is 



2t 


.a 












n 




( ■> “g 
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“ '’g 1 


1 r dr dS 





2irr 




Q + 
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M 



M 



Gr 



Q 



Gr 



w 



(4-2-1) 

b 



which reduces to 



a 

I q Wq - w j r dr = 
b 



( -w 



a 



G ^ 



*G ^ - 



M- i 
Gr 



Qf, w 

Gr 



(4-2-2) 



So far, it has only been asserted that the integral 
equation for the circular plate case should have a form similar 
to Eq, (4-2-2). To obtain the exact integral equation, a 
direct method is used. 



As a starting point for the direct method, consider the 
following expression: 



U 





dr 



This expression can be written slightly differently as V 
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3 J ^ J ^ 
d w d w. 



2 2 
dr dr 



r dr 



Integrating both expressions by parts yields 
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dr 
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Recalling Eq. (2-3-2), it is possible to rewrite these 
expressions in terms of shear Qj. 
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Setting U equal to V and canceling some terms gives 
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dw^ d w 
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The two integral terms above are integrated by parts, with 
the operator of Eq. (2-3-5) used to simplify the expressions 



a 


dw Q 
G r 

— — r 


dr 


_ 


Q 

r 

w r - 


a 


b 


dr D 






D 


b 


a 


dw Q 

r 


dr 




®Gr 
w r — 


a 


b 


dr D 






D 


b 



a 



V w w^ r dr 
G 



V w_ w r dr 
G 



which, when substituted into the previous equation, yield 
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After some simplification, this expression becomes 
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The last equation bears a close resemblance to Eq. (4-2- 
2), except for sign differences and the fact that terms 
involving the second derivative of deflection have yet to be 
replaced by terms involving moment. To accomplish the latter, 
the following expression is added to the third term of the 
right hand side of the preceding expression, and subtracted 
from the fourth term 

V dw dw„ 

^ Cj 

r dr dr 

Through Eq. (2-1-7) this gives 
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(4-2-3) 



which can be rewritten using Eq. (2-3-4) as 
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(4-2-4) 



Remarkably, with the exception of a sign difference that is 
attributable to the difference in sign convention, this 
expression is identical to Eq. (4-2-2), which was asserted 
based on the integral equation developed by Butterfield. 
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4.3 Selection of a Greenes Function 



As is the case for the beam described in Chapter 3, Eq. 
(4-2-4) is useful only if wq is chosen as a Green’s function 
with certain properties. As will become evident later in this 
chapter, the Green’s function must describe a ring loaded plate 
to solve circular plate bending problems, just as a Green’s 
function describing a point loaded beam was used to solve the 
beam bending problem. 

To solve the beam problem, Butterfield choses a Green’s 
function describing an infinitely long beam with a point load. 
It is not essential that a Green’s function for an infinite 
beam be used. In fact, as one would expect from Betti’s 
reciprocal work theorem, a Green’s function describing a finite 
geometry will work as well, provided the function also 
describes the response to a point load. 

Applying Butterfield’s method to the circular plate case, 
a Green’s function representing the ring loaded plate shown in 
Figure 4-2 is adopted. The Green’s function and associated 
derivatives are included in Appendix A. It should be noted 
that the outer and inner radii of the Green’s function (c and 
d) need not coincide with the outer and inner radii of the 
actual plate (a and b). 
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Figure 4-2. Plate geometry for the selected Green’s 
function 



That a simple support is chosen for the outer plate radius is 
not important; a plate clamped at both the inner and outer 
radii would also have been suitable. What is important is that 
the loading condition for the Green’s function case be a ring 
load. To simplify calculations, a ring load of magnitude 1 is 
used . 



4 . 4 Solution of Integral Equation 

Having selected a Green’s function, the next step is to 
obtain a form of Eq. (4-2-4) suitable to solve the actual plate 
bending problem. Rearranging terms, Eq. (4-2-4) can be 
rewritten as . 
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The ring load may be expressed in terms of the dirac delta 
function 6(r,rQ) 

Substituting this expression into Eq. (4-4-1) and making use of 
the following property of a dirac delta function 
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(4-4-2) 



By choosing Cq at radius a, and then at radius b, the 
deflection w in the left hand side of Eq. (4-4-2) is the 
deflection at the boundary, which is either a known or unknown 
boundary condition. Hence, two equations are available to 
solve for the four unknown boundary conditions. 
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To obtain two more equations, Eq. (4-4-2) is 
differentiated with respect to r^. This yields the second 
Green’s function, which is evaluated with r^ at a and then 
again at b. The second Green’s function corresponds to the 
slope of the plate, which at r=a and r=b is again either a 
known or an unknown boundary condition. There is thus now a 
complete set of four equations to solve for the four unknown 
boundary conditions. 



Eq. (4-4-2) can be rewritten as 
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(4-4-3) 



Differentiating this expression with respect to ro gives 
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Eqs. (4-4-3) and (4-4-4), evaluated both for ro=a and ro=b, 
give rise to a system of four equations that when solved 
simultaneously, yield the four unknown boundary conditions, 
matrix form, these equations can be written as 
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Using Gauss elimination, the preceding system of equations is 
solved for the four unknown boundary conditions. 



The final step is to calculate deflection, slope, moment 
and shear across the plate. Deflection and slope are found 
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using Eqs. (4-4-3) and (4-4-4) by substituting in the 
calculated boundary conditions and evaluating the expressions 
at the value of ro of interest. Moment and shear are 
determined using Eqs. (2-1-7) and (2-3-2), which require 
calculation of the second and third derivatives or w with 
respect to ro* Differentiating Eq. (4-4-4) results in the 
desired expressions 
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making a complete solution for the circular plate bending 
problem in the elastic range possible. 
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CHAPTER 5 



NON-LINEAR BEHAVIOR 



This chapter develops a boundary integral formulation for 
circular plate bending problems that includes non-linear 
material behavior. Specifically, the boundary integral 
equation developed in Chapter 4 for the elastic plate is 
modified to include plastic moment terms to account for the 
plastic behavior. An incremental load method similar to that 
used by Moshaiov and Vorus can then be applied to arrive at a 
solution to the plate bending problem once yielding has begun. 

5 . 1 Plastic Stress-Strain Relationships 

When considering plasticity, strain can be thought of as 
having two components: an elastic component and a plastic 

component. In cylindrical coordinates, which is applicable to 
the two-dimensional circular plate case, the total strain e can 
be expressed in terms of these two components as 



e ^ p 

€ = £ + £ 
r r r 



^ P 

"e - ^0 ^ ^0 



(5-1-1) 



where the superscripts e and p refer to the elastic and plastic 
strain components, respectively. 

Eq. (2-1-4), which provides the relationship between 
stress and strain for a linearly elastic material, can thus be 
rewritten in terms of the total strain and plastic strain as 
follows : 
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5 . 2 Plastic Moment-Curvature Relationships 

To develop the moment-curvature relationships with 
plasticity, Eq. (5-1-2) is substituted into Eq. (2-1-3) to give 
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These expressions can be rewritten as 
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where the radial and tangential plastic moments, Mj.P and MgP, 
are defined as 
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Recalling the relationships for total strain of Eq. (2-1-6) and 
performing the integration over the plate thickness yield the 
moment-curvature relationship for a circular plate with 



47 



plast icity 
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5 . 3 Governing Differential Equation With Plasticity 

Chapter 2 develops the governing differential equation for 
a circular plate in the elastic range by substituting the 
elastic moment-curvature relationships into the elastic 
equilibrium equation. In this section, a similar approach is 
used to develop the governing differential equation with 
plasticity. The difference here is that the moment-curvature 
relationships (Eq. (5-2-3)) include plastic terms. 

Substituting Eq. (5-2-3) into the equilibrium equation for 



a circular 


plate 


(Eq. 


(2- 


3-1 ) ) gives 
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= (5-3-1) 

D 

where the tilda is used to distinguish the fact that the shear 
includes plastic effects. Differentiating this expression once 
with respect to r gives the fourth order governing differential 
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equation for the plate with plasticity 



4 3 

d w 2 d w 

dr ^ r dr ^ 



2 

1 d w 1 dw 

Z ^ ^ 7 ~ 

r dr r dr 



q 

D 



1 

D 



2 dM P d^M P 

^ ^ 

r dr dr^ 



1 dM^ 
r dr 



which can be rewritten using the operator defined in Eq. (2-3- 
5 ) as 
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This equation is identical to the equilibrium equation for the 
elastic case with the exception of the terms that include 
derivatives of the plastic moments. Note that these terms may 
be thought of as additional pseudo loads that must be combined 
with the actual load to account for the plastic behavior. 



5 . 4 Boundary Integral Formulation With Plasticity 

In the preceding chapter, the governing differential 
equation is used to develop a boundary integral formulation for 
a circular plate in the elastic case. This section presents a 
similar approach to develop a boundary integral formulation 
for the plastic case. 

As a starting point, Eqs. (2-3-4) and (5-3-2) are 
substituted into Eq. (4-2-3) to give 
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This equation is not yet in a form that is usable for 
solving the circular plate bending problem. In the first 
place, there are several terms that include derivatives of the 
radial or tangential plastic moments, which are not generally 
known across the plate. Secondly, the shear and moment terms 
Qj. and Mj- evaluated at the boundaries are not the actual shears 
and moments, since they were developed based on elastic 
cons idel*at ions only. 

To resolve these concerns, the integral term involving 
plastic moments in Eq. (5-4-1) is integrated by parts as 
follows : 



50 





dM P 


d^M P 


dM-P 




2 — ^ + 


r 

r 


0 


Vj 

h 


dr 


dr^ 


dr 

j 



dr 



M, 



dM ^ 

M P + r — - M P 
dr ® 



dw. 



dr 



dM * 

M P + r — 

^ dr 



- M. 



dr 



Integrating by parts a second time and collecting terms gives 
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Using the above relationship and choosing a Green’s function 
that describes a ring loaded plate as is done for Eq. (4-4-2), 
Eq. (5-4-1) becomes 
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Recall that and Mj. represent shear and moment for the 
elastic case. However, they are not the actual shear and 
moment for the plastic case. With plasticity, the actual 
shears and moments are combinations of elastic and plastic 
terms. From Eq. (5-3-1), the shear relationship with 
plasticity is 
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and from Eq. (5-2-3), the moment relationship with plasticity 
is 
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where, as with the case for shear, the tilda is used to 
distinguish the fact that the moment term includes plastic 
effects. When substituted into the preceding expression, these 
expressions give the desired boundary integral formulation for 
the circular plate bending problem 
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By evaluating Eq. (5-4-2) with the ring load positioned 
both at ro=a and ro=b, two equations are obtained to solve for 
the unknown boundary conditions. Eq. (5-4-2) is differentiated 
with respect to r<j to obtain a second Green’s function. This, 
via Eq. (2-1-1), gives an expression for slope of the actual 
plate. The second Green’s function is evaluated with the ring 
load positioned at r <3 = a and ro = b, to provide the remaining two 
equations . 



Differentiating Eq. (5-4-2) with respect to rg gives 
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In matrix form, the system of equations to be solved for 
the four unknown boundary conditions for the non-linear case 
can be written using the notation of Eq. (4-4-5) as follows; 
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Observe that Eq. (5-4-4) is identical to Eq. (4-4-5) 
except for the additional moment terms in the last matrix of 
Eq. (5-4-4). These additional terms account for the plastic 
behavior. 



To find moments and shears across the plate, expressions 
for the second and third derivatives of deflection with respect 
to ro are required. These expressions follow: 
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CHAPTER 6 



NUMERICAL SOLUTION 



This chapter describes the numerical procedure that was 
used to solve the circular plate bending problem in both the 
elastic and plastic ranges. The incremental load approach used 
in the plastic range is first described. A brief description 
of the computer program that was developed based on the 
equations of Chapters 4 and 5 follows, including a discussion 
of some aspects of numerical methods used in this program. 

6 . 1 Incremental Load Method 

Once plastic deformation starts to occur in the plate, the 
linear relation between stress and strain is not valid. As 
discussed by Mendelson^^^ , strains in the plastic range are no 
longer uniquely determined by the stresses, and in fact depend 
on the loading history. Therefore, it is not possible to use 
the relationships developed in Chapter 5 to directly arrive at 
a plate bending solution given a loading condition outside the 
elastic range. Rather, an incremental approach must be used 
wherein the load is increased in small steps once yielding 
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occurs at any point in the plate, and the complete stress state 
for the plate is determined before load again is increased. 

The incremental load method used for this analysis is 
based on that adopted by Moshaiov and Vorus, and is frequently 
seen in finite element solutions of structural analysis 
problems. The principal steps of this method are outlined 
below: 

1. Using the plate parameters and loading conditions for 
the bending problem at hand, an elastic solution is obtained 
with the equations developed in Chapter 4. 

2. The load at which the onset of yielding will occur in 
the plate is calculated based on the stress state determined in 
Step 1. The selected yield criterion is the von Mises 
criterion, which for the symmetrically loaded circular plate 
case has the form: 




where Oq is the equivalent, or effective, stress which 
represents the von Mises yield surface. Yielding will occur 
when Oq equals or exceeds the uniaxial yield stress (Sy). 

3. Unknown boundary conditions at yielding are 
determined using Eq. (5-4-4), with plastic moments initially 
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set equal to zero. Total strains at predetermined integration 
points are calculated using Eqs. (2-1-6), (5-4-3) and (5-4-5), 

From Eq. (2-1-4), the stress state at the integration points at 
yielding is determined, and is stored to be later updated as 
load is increased. Also stored for later updating are the 
deflections, slopes, moments and shears across the plate. 

4. The load which caused yield is then increased by a 
small incremental amount. The increments of the total strains 
resulting from this incremental load are calculated using Eq. 
(2-1-6) in an incremental form, from which elastic strains can 
be calculated by means of Eq. (5-1-1). 

5. Using Eq. (2-1-4), the elastic stress increment 
corresponding to the applied incremental load is calculated. 
This incremental stress is added to the stress stored in Step 3 
to give the total stress at the end of the load increment. A 
check is made throughout the plate using the yield criterion of 
Step 2 to determine where yielding has occurred. 

6. In those plate regions where yielding has occurred, 
the plastic strain increment (that is, the plastic strain due 
to the incremental load) is calculated. Only materials that 
exhibit e las t ic-perf ect ly plastic behavior as depicted in 
Figure 6-1 have been examined. However, the method can be used 
with other material behavior, such as strain hardening. For 
materials exhibiting elastic-perfectly plastic behavior, the 
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plastic strain increment 6cp is the total incremental strain 
determined from Eq. (2-1-6). This is shown in Figure 6-1. 

7. The plastic moments are determined using Eq. (5-2-2). 

8. Steps 3 through 6 are repeated, with the plastic 
moments calculated in Step 7 used as an input to Eq. (5-4-4). 
Step 7 is repeated and the plastic moments compared to the 
moments calculated previously in Step 7. Iterations are 
performed as necessary until these values converge, according 
to a predetermined convergence criteria. 

9. After convergence of plastic moments is achieved in 
Step 8, stresses, plastic strains, deflections, slopes, moments 
and shears throughout the plate are updated. 

10. Another increment of load is applied, and Steps 3 
through 9 above are repeated, using the plastic moment from the 
previous load step as an initial estimate for plastic moment. 
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Figure 6-1. Uniaxial stress-strain curve for eiastic- 
perfectly plas-tic material. 



6 . 2 Computer Program Description 

Computer programs written in the Fortran 77 language have 
been developed to solve the circular plate bending problem and 
are included as Appendix B. Three programs, supported by nine 
subroutines, are used. 

The programs INPUT and LOAD create data files which are 
used by the program MAIN to solve the bending problem. INPUT 
and LOAD prompt the user for information describing the problem 
to be solved, such as material properties, plate geometry, 
boundary conditions and plate loading conditions. In addition, 
certain adjustable parameters are input, which include step 
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sizes for the load increments, number of integration 
increments, and the desired percentage for convergence of the 
plastic moment increments. The program MAIN determines the 
plate bending solution, both in the elastic and plastic ranges. 

A flow diagram outlining the major elements and overall 
logic path of this program is shown in Figure 6-2. 
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6 . 3 Numerical Methods 



Eqs. (5-2-1), (5-4-4), (5-4-5) and (5-4-6) require 

integrations to be performed across the plate. These 
integrations are accomplished by a trapozoidal rule numerical 
integration scheme. It was found that 10 increments across 
both the plate radius and the plate half thickness gave 
satisfactory results. 

To ensure convergence of plastic moments, a root mean 
square average of the the plastic moments at each station 
across the plate is calculated and compared to the root mean 
square average from the previous iteration. An agreement of 
less than 0.1% was used for the analysis work for which results 
appear in Chapter 7. 
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CHAPTER 7 



RESULTS AND DISCUSSION 



This chapter presents results of analyses performed on 
three different plate configurations using the computer program 
described in the preceding chapter. Results in the elastic 
range for all three cases are compared with the analytic 
solution from Roark^^^, and for one case in the plastic range 
where published results using another solution method are 
available. Descriptions of the case studies are given in 
Sections 7.1 through 7.3. A discussion of the results is given 
in Section 7.4. 

7 . 1 Simply Supported Elasto-Plastic Annular Plate 

The annular plate shown in Figure 7-1 is first examined. 
The plate is simply supported at both the inner and outer 
radii, and is subjected to a uniform lateral load. A material 
exhibiting elast ic-perf ect ly plastic behavior as depicted in 
Figure 6-1 is used. The yield stress (Sy) is 16 ksi and the 
Young’s modulus is 10x10^ ksi. Material properties and plate 
dimensions were selected to allow comparison with results 



66 



obtained using a finite element method by Armen et 

Plate deflection at yielding is shown in Figure 7-2, which 
also includes results from the analytic solution given by Roark 
for comparison. Plate deflections for the present solution at 
several loads in the plastic range are shown in Figure 7-3, as 
well as results obtained by Armen et al using a finite element 
method of analysis. 

Figure 7-4 shows the plastic zones in the plate at the 
three load conditions depicted in Figure 7-3. Results obtained 
by Armen et al are also included. 



Sy= 16 KSI 

E = 10X10^ KSI 





Figure 7-1. Simply supported annular plate problem 
description . 
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,ure 7-2. Deflection of simply supported annular pLa at 




RADIAL POSITION IN INCHES 

Figure 7-3. Deflection of simply supported annular plat in 
the plastic range. 
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"igure 7-4. Plastic zones in simply supported annular pla"j. 

7 . 2 Clamped Elasto-Plast ic Annular Plate 

The annular plate shown in Figure 7-5 is next examined. 

The plate is clamped at both the inner and outer radii, and is 
again subjected to a uniform lateral load. A material 
exhibiting the same properties as described in Section 7.1 is 
used. For the clamped annular plate case, results obtained 
using other methods could not be found for the plastic range. 
Consequently, material properties and plate dimensions were 
selected to allow a qualitative comparison with results for the 
simply supported plate of Section 7.1 

Plate deflections at the onset of yielding are shown in 
Figure 7-6, and at several loads in the plastic range in Figure 
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7 7. The deflection at yielding using the analytic solution 
given by Roark is shown in Figure 7-6 for comparison. Finally, 
Figure 7-8 shows the plastic zones in the plate at the three 
load conditions depicted in Figure 7-7. 



Sy= 16 KSI 
E = 10X10^ KSI 
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Figure 7-5. Clamped annular plate problem descript icn . 
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re 7-6. Deflection of clamped annular plate at the 
onset of yielding 




Figure 7-7. Deflection of clamped annular plate in the 
plastic range 
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Figure 7-8. Plastic zones in clamped annular plate. 

7 . 3 Simply Supported/Guided Elasto-Plastic Annular Plate 

The third annular plate to be examined is shown in Figure 
7-8. The plate is simply supported at the outer radius, with 
the inner radius guided. The plate is again subject to a 
uniform lateral load. A material exhibiting the same 
properties described in Section 7-1 is assumed. 

For this case, results obtained using other methods could 
again not be found for the plastic range. However, a solution 
for a continuous simply supported circular plate has been 
obtained by Moshaiov and Vorus, and offers a qualitative 
comparison with the results for the simply supported/guided 
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annular plate. 



Plate deflection at the onset of yielding is shown in 
Figure 7-10 and at several loads in the plastic range in Figure 
7-11 for the present analysis method. The deflection at 
yielding obtained using the analytic solution from Roark is 
shown in Figure 7-10 for comparison. Figure 7-12 shows the 
plastic zones in the plate at the three load conditions 
depicted in Figure 7-11 for the present solution. Corres- 
ponding results from Moshaiov and Vorus for the continuous 
plate of Figure 7-13 are shown in Figures 7-14 and 7-15. 



Sy= 16 KSl 
E = 10X10^ KSI 




Figure 7-9. Simply supported/guided annular plate probl^.u 
description . 
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radial position in inches 



Figure 7-10. Deflection of simply supported/guided annular- 
plate at the onset of yielding. 




radial position in inches 



Figure 7-11. Deflection of simply supported/guided annular 
plate in the plastic range 
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7 . 4 Discussion of Results 



The results presented in Sections 7.1 through 7.3 for the 
elastic range are uniformly in excellent agreement with the 
results using the analytic solutions given by Roark. This is 
as should be expected, since the elastic range solution using 
the present method is based on a closed form analytic solution. 

Results in the plastic range, while reasonable, are not in 
as close agreement with other published results. The only case 
where data for a direct comparison could be found is the 
configuration discussed in Section 7.1. From Figure 7-3, the 
deflections at loads of 652 psi and 852 psi for the present 
solution are in excellent agreement with the finite element 
results of Armen et al. 

Results for the annular plate of Section 7.1 at the 1052 
psi load are not as encouraging. Deflections and the extent of 
the plastic zone are greater for the solution obtained by Armen 
than for the present solution. To explain this difference, the 
size of the integration and load increments has been varied, 
but it shows little effect on the results. The reason behind 
this difference is not known. However, it is believed that it 
might be due to a programming error either in this report or in 
the reference. Also, differences could arise by not taking 
into account additional integration terms arising from the 
singularity at r=r^, as discussed in Appendix A. 
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Finally, it is of interest to compare the results for the 
simply supported/guided plate of Figure 7-9 to the results for 
a simply-supported circular plate obtained by Moshaiov and 
Vorus. It is recognized that the two cases are not equivalent. 
However, one would expect the two plates to behave globally in 
a similar fashion. This is in fact the case, with both 
deflections and plastic zones showing reasonably good agreement 
in view of the differences between the two configurations. 



CHAPTER 8 



SUMMARY AND CONCLUSIONS 



8 . 1 Summary 

This thesis has developed a formulation for solving 
axisymmetrically loaded annular plate bending problems using 
boundary integrals. This particular formulation is unique in 
that it treats the annular plate as a one-dimensional problem, 
using "ring" type Green’s functions to determine unknown 
boundary conditions and arrive at a plate bending solution. 

The formulation allows a numerical incremental load method to 
be used in the plastic range for the treatment of non-linear 
behavior. Results in the elastic range show excellent 
agreement with results obtained using a conventional analytic 
solution. In the plastic range, results are in reasonable 
agreement with those obtained using a finite element method. 

8 . 2 Conclusions 

This thesis treats the two-dimensional problem of analysis 
of axisymmetrically loaded annular plates as a one-dimensional 
problem. Using a method of solving one-dimensional problems 
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with boundary integrals developed by Butterfield, the thesis 
demonstrates that a closed form solution for the annular plate 
problem can be obtained. It further demonstrates that this 
approach is suitable for the use of an incremental load method 
to obtain a solution in the plastic range. 

The formulation developed in this thesis is advantageous 
over a conventional analytic solution in that it allows a wide 
range of problems with different loading and boundary 
conditions to be solved using a single algorithm. Further, the 
simplicity of the approach is useful from an educational 
standpoint in the teaching of boundary element methods. 

8 . 3 Recommendations for Future Work 

There remains room for much additional work in this area. 
Some suggestions follow: 

1. A formulation for axisymmetrically loaded continuous 
circular plates remains to be developed. The work presented 
here lays a foundation for the continuous plate formulation. 
However, extending this method to the continuous plate case may 
require some modifications. 

2. In developing a formulation for the simple beam 
problem using boundary integrals, Butterfield uses both a 
direct and an indirect method, and shows that they yield the 
same result. Only the direct method is presented in this 
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thesis. The indirect method is a more intuitive approach than 
the direct method, and offers advantages in furthering the 
understanding of a boundary integral solution. It would, 
therefore, be worthwhile to pursue an indirect method as well 
as a direct method in developing a boundary integral 
formulation for circular plate geometries. 

3. The Green’s functions selected for this anlysis 
involve many terms and are awkward from a computational 
standpoint. More thought needs to be given to selecting a 
"ring" type Green’s function that is simpler and therefore 
offers advantages in simplifying calculations and reducing 
computational time. 

4. For the plastic range, the lack of close agreement 
between the results presented in this work and those obtained 
by Armen et al using a finite elemet method need to be better 
understood. Additional comparative data should be found or 
developed to increase confidence in the results presented in 



this work. 



APPENDIX A 



SELECTED GREEN’S FUNCTIONS 



This Appendix presents the selected Green’s functions and 
required derivatives that are used to solve the plate bending 
problems in this analysis. The functions which mathematically 
describe the deflection, slope, radial and tangential moments, 
and shear are obtained from Roark. 

A . 1 Description of Geometry and Sign Conventions 

As discussed in Chapter 4, the first Green’s function 
selected for this analysis describes the response of annular 
plate that is simply supported at the outer radius and 
unsupported along the inner radius to a ring load. Figure A-1 
depicts this plate configuration and important notation. 
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Figure A-1. Representation of Green’s function used in 
the analysis . 



Roark uses a sign convention with deflection wq positive 
upward, slope iQ positive when the deflection wq increases 
positively as r increases, moment Mj. positive when creating 
compression on the top surface of the plate, and the shear 
force Qgt positive when acting upward on the inner edge of an 
annular section. Subscripts c and d refer to the radial 
pos it ion . 

This sign convention differs from the sign convention of 
Timoshenko, which has been adopted for this analysis. 
Accordingly, adjustments must be made. Specifically, 
Timoshenko takes deflection to be positive downward and shear 
to be positive when acting downward on the inner edge of the 
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plate, such that the signs for these items as given by Roark 
must be changed. 

For slope, Roark uses a convention opposite to that of Eq. 
(2-2-1) and the convention used by Timoshenko. This 
difference, combined with the difference in defining the sign 
of deflection, effectively cancel each other, so the sign of 
slope from Roark does not change for this analysis. Finally, 
for moments, Roark and Timoshenko use the same sign convention, 
so no sign changes for moment terms are necessary. 

The formulas in the sections which follow include a number 
of general plate functions and constants that are not depicted 
in Figure A-1. Functions L and G and their derivatives are 
included in Section A. 7; constants F and C are included in 
Sect ion A. 8 . 

A . 2 Singularities 

It is not possible to evaluate the Green’s function and 
associated derivatives at the exact radial location where the 
ring load is applied, due to singularities. To avoid this 
problem on the boundaries, the ring load is positioned a small 
distance inside the outer plate radius and outside the inner 
plate radius (0.00001 inches for the results presented in 
Chapter 7). During final review of this thesis, it was noted 
that because of the singularity at r = r< 5 , the domain integrand 
involving plastic moments of Eqs. (5-4-2) through (5-4-6) 
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includes singularities. Therefore, the limiting values should 
be explored and a correction applied in way of the singular- 
ities. Insufficient time was available to repeat the analysis 
with the applied correction to see its effect on the results. 
Note that this correction affects only the solution in the 
plastic range. 



A . 3 Deflection. Slope. Moment and Shear 

The following expressions describe deflection wq , slope 
radial moment Mqj.* and shear Qq corresponding to the first 
Green’s function: 
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The expression <r-rQ>° represents a singularity function. 
The function is equal to zero if r<rQ. If r>rQ, the brackets 
become like any other brackets. Hence, for r>ro» <r-ro>° is 
(r-rg) to the power of 0, which equals to 1. 



A . 4 First Derivative of Deflection. Slope, Moment and Shear 
The expressions which follow describe the first deriv- 
atives with respect to rg of wq, Mqj. and QQr* first 

derivative of wq corresponds to the second Green’s function 
required by the analysis. 
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A . 5 Second Derivative of Deflection, Slope, Moment and Shear 
The following expressions describe the second derivatives 
with respect to r^ of wq , Mq^ and Qq^ 
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A . 6 Third Derivative of Deflection, Slope, Moment and Shear 
The following expressions describe the third derivatives 
with respect to ro of wq , and Qq^ 
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A . 7 Plate Functions L3, L9, G3> G6 and G9 and Their 

Derivatives 
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A . 8 Plate Constants Fi. F 4 . F?. Ci. and C 7 
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A . 9 Derivatives of the Greenes Function With Respect to r 
The equations describing non-linear plate behavior 
included in Section 5-4 are in terms of the first and second 
derivatives of the Green’s function wq with respect to r. 

These derivatives can be expressed with the aid of Eqs. (2-1-1) 
and (2-1-7) in terms of moments and slopes as follows: 
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Expressions for slopes and moments for the Green’s 
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function (iq and Mqj.) are available and have been presented in 
Section A. 3. It is further noted that Eqs. (5-4-2) through (5- 
4-6) require derivatives of the above two expressions with 
respect to rQ. This is accomplished via the expressions for 
derivatives of ♦g and Mqj. with respect to ro included in 
Sections A. 5 through A. 7 above. 



93 



APPENDIX B 



COMPUTER PROGRAMS 



This Appendix contains the computer programs developed to 
perform the analysis work. A flow chart showing the overall 
program structure is presented in Chapter 6. 

The Appendix is organized with the programs INPUT and LOAD 
presented first, followed by the program MAIN. The supporting 
subroutines for the program MAIN follow that program. 
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c 

c* 

c* 

Ct PROGRAM INPUT 
C* 

C* THIS PROGRAM ALLOWS THE USER TO DEFINE THE PROBLEM OF INTEREST 
C» AND TO SET CERTAIN ADJUSTABLE PARAMETERS SUCH AS NUMBER OF DEPTH 
C* INTEGRATION INCREMENTS. LOADING CONDITIONS ARE INPUT USING THE 
C» PROGRAM LOAD. TO RUN, THIS PROGRAM MUST ACCESS THE FILES GEOM.DAT, 
C* BC. DAT, BCPOS.DAT AND ADJ.DAT 
Ct 

Cttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

c 

PROGRAM INPUT 
C 

IMPLICIT REALMS (A-H,0~Z) 

CHARACTER*! BFLAGl (4) , BFLAG2(4) , RKSP 
CHARACTER*19 LBL( 8) , LGEOM(9 ) 

CHARACTER*25 LADJ(6) 

C 

RBAL*8 ADJ(7) ,BC(4) ,GBOM(9) ,NU 
C 

INTEGER POS(8) ,BCP0S(8) ,M,N 
C 

C* 

C****************************************************************** 

c* 

C* INPUT MATERIAL CONSTANTS AND PLATE GEOMETRY. THE FLEXURAL 
C* RIGIDITY IS CALCULATED BASED ON THE INPUT PARAMETERS. N IS THE 
C* TOTAL NUMBER OF PARAMETERS 
C» 

c****************************** ************************************ 

c* 

c 

N=9 

C 

LGBOM(l)-*OUTER RADIUS A = * 

LGEOM(2)=* INNER RADIUS B = * 

LGB0M(3)=’P0ISS0N RATIO NU = * 

LGB0M(4)=*PLATB CONSTANT D = * 

LGBOM(5)=*THICKNBSS T = * 

LGEOM(6)=* YIELD STRESS = * 

LGEOM(7)=* YIELD STRAIN = ’ 

LGEOM(8)= ’ ULTIMATE STRESS = ’ 

LGEOM(9)=*ULTIMATE STRAIN = ’ 

C 

WRITE(*, 1) 

1 FORMAT(///,6X, ’WOULD YOU LIKE TO SEE THE CURRENT PLATE DIMENSIONS 
,/,2X,’AND MATERIAL CONSTANTS (Y/N) ? *\) 

READ(*,30) RESP 
C 

IF (RESP .EO. ’Y*) THEN 
WRITE(*,2) 
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QO O) 



2 FORMAT(//,2X, ’THE CDRRBNT VALDES FOLLOW. NOTE THAT A ZERO’, 

+ /,2X, ’VALUE IS SHOWN FOR D WHEN T AND E ARE GIVEN AND VICE’, 

+ /,2X, ’ VERSA ...’,/) 

OPEN(6, FILE= ’GEOM.DAT’ , STATUS =’ OLD ’ ) 

WRITE(*,30) ’ ’ 

DO 8 1=1, N 

READ(6,70) GBOM(I) 

WRITB(*,6) LGBOM(I),GBOM(I) 

FORMAT(2X, A22,F20. 10) 

CONTINUE 
CL0SE(6) 

ENDIF 
C 

WRITE(*,9) 

9 FORMAT(//,6X, ’DO YOU WISH TO INPUT OR CHANGE GEOMETRY OR*, 

+ /,2X, ’MATERIAL CONSTANTS (Y/N) ? ’,\) 

READ(*,30) RESP 
C 

IF (RESP .EO. ’Y’) THEN 
C 

WRITE(», 10) 

10 FORMAT(/,2X, ’INPUT THE VALUES USING DECIMAL NOTATION ’,/) 

DO 12 1=1,3 

WRITE(»,20) LGEOM(I) 

RBAD(»,70) GEOM(I) 

12 CONTINUE 
C 

DO 13 1=5,9 

WRITB(», 20) LGEOM(I) 

READ(*,70) GEOM(I) 

13 CONTINUE 

GBOM(4)=GEOM(6)/GEOM(7)»GEOM(5)**3. 0/12. 0/( 1 .0-GEOM(3)»*2. 0) 

C 

OPEN(6,FILE=* GEOM.DAT’ , STATUS= ’ NEW ’ ) 

DO 15 1=1, N 

WRITE(6,70) GEOM(I) 

15 CONTINUE 

CLOSE(6) 

ENDIF 

C 

C 

Ctt*****t******ttt**t************tt**********t****t*************t** 

c 

C INPUT PLATE BOUNDARY CONDITIONS 
C 

Cttttttttt*tt***t**t*tttt*t*ttt*********t*t****ttt*******t**tttt**t 

c 

c 

LBL(1)=’SHEAR AT A = ’ 

LBL(2)=’MOMENT AT A = ’ 

LBL(3)=’SLOPE AT A = ’ 

LBL(4)=’DEFLECTION AT A = ’ 

LBL(5)=’SHEAR AT B = ’ 

LBL(6)=’MOMENT AT B = ’ 

LBL(7)=’SLOPE AT B = ’ 
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c 



LBL(8)=*DEFLBCTI0N AT B 



0PEN(14,PILE=’BCLBL. DAT ’ , STATUS= ’ NEW’ ) 

DO 155 1=1,8 

WBITE(14, 156) LBL(I) 

155 CONTINUE 

156 POBMAT(A22) 

CLOSE( 14) 

C 

WBITE(»,160) 

160 F0RMAT(//,6X, 'WOULD YOU LIKE TO SEE THE CUBBENT PLATE BOUNDABY’ 

+ ,/,2X, 'CONDITIONS (Y/N) ? ’\) 

BEAD(*,30) BESP 
C 

IF (BESP .EQ. 'Y') THEN 

OPEN(ll,FILE=’BC. DAT’ , STATUS= ’ OLD ’ ) 

OPEN ( 12, FILE=' BCPOS.DAT* , STATUS= * OLD ’ ) 

WBITE(»,30) ’ * 

DO 165 1=1,4 

HEAD(11,70) BC(I) 

165 CONTINUE 

DO 190 1=1,8 

BEAD(12,100) BCPOS(I) 

IF (BCPOS(I) .GT. 0) THEN 

WBITE(*,186) LBL(I) ,BC(BCPOS(D) 

186 FOBMAT(2X,A22,F20. 10) 

ENDIF 

190 CONTINUE 

CLOSE( 11) 

CLOSE( 12) 

ENDIF 

C 

WHITE(», 17) 

17 F0HMAT(//,6X, 'DO YOU WISH TO INPUT OB CHANGE PLATE BOUNDARY’, 

+ /,2X, ’CONDITIONS (Y/N) ? ’,\) 

READ(»,30) BESP 
C 

IF (BESP .EQ. ’Y*) THEN 
WBITE(», 19) 

19 F0RMAT(//,6X, 'INDICATE WHETHER THE BOUNDARY CONDITIONS FOR THE’ 

e /,2X, ’PLATE ARE KNOWN OR UNKNOWN BY TYPING IN THE LETTER "U"’ 
e /,2X,’FOR UNKNOWN BOUNDARY CONDITIONS AND HITTING RETURN FOR’ 
e /,2X, ’KNOWN BOUNDARY CONDITIONS’,//) 

C 

WRITE(«,20) 'SHEAR AT A ? ’ 

READ(*,30) BFLAGl(l) 

WRITE(»,20) 'MOMENT AT A ? ’ 

READ(*,30) BFLAGK2) 

WRITE(»,20) 'SLOPE AT A ? ’ 

BEAD(A,30) BFLAGK3) 

WRITE(»,20) 'DEFLECTION AT A ? ’ 

BEAD(*,30) BFLAGK4) 

C 

WRITE(»,20) 'SHEAR AT B ? ’ 

READ(*,30) BFLAG2(1) 
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WHITE(*,20) ’MOMENT AT B ? ’ 

READ(>,30) BFLAG2(2) 

WRITE(»,20) ’SLOPE AT B ? ’ 

READ(*,30) BFLAG2(3) 

WRITE(>,20) ’DEFLECTION AT B ? ’ 

READ(*,30) BFLAG2(4) 

20 FORMAT(2X»A22,\) 

30 FORMAT(Al) 

C 

WRITE(»,40) 

40 FORMAT(//,6X, ’TYPE IN THE KNOWN BOUNDARY CONDITIONS USING 
«2X. ’DECIMAL NOTATION WHEN PROMPTED’,//) 

C 

OPEN(8, FILE=’ BFLAG1.DAT’ ,STATUS= ’NEW’ ) 

OPEN( 9, FILE=’ BFLAG2.DAT’ , STATUS =’ NEW ’ ) 

OPENdO, FILE=’BC.DAT’ , STATUS^ ’ NEW ’ ) 

C 

DO 50 1=1,4 

WRITE(8,30) BFLAGl(I) 

IF (BFLAGl(I) .NE. ’U’) THEN 
L=L+1 

WRITE(*,20) LBL(I) 

READ(»,70) BC(L) 

WRITE(10,70) BC(L) 

ENDIF 

50 CONTINUE 
C 

DO 60 1=1,4 

WRITE(9,30) BFLAG2(I) 

IF (BFLAG2(I) .NE. ’U’) THEN 
L = L + 1 

WRITE(>,20) LBL(I+4) 

READ(>,70) BC(L) 

WRITE(10,70) BC(L) 

ENDIF 

60 CONTINUE 
C 

70 FORMAT(F20. 10) 

C 

CLOSE(8) 

CLOSE(9) 

CLOSE( 10) 

C 

OPEN(ll , FILE=’POS.DAT’ , STATUS= ’ NEW ’ ) 

OPEN( 12, FI LE=’ BCPOS.DAT’ , STATUS =’ NEW ’ ) 

C 

L = 0 

DO 80 1=1,4 

IF (BFLAGKI) .EQ. ’U’) THEN 
L = L+1 
POS(I)=L 
BCPOS( I) =0 
ELSE 
M=M+1 

BCPOS( I)=M 
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POS( I)=0 
ENDIF 
C 

WRITEdl, 100) POS(I) 

KRITE(12. 100) BCPOS(I) 

C 

80 CONTINUE 
C 

DO 90 1=5,8 

IF (BFLAG2(I-4) .EO. *U») THEN 
L=L-^1 
POS(I)=L 
BCPOS(I)=0 
ELSE 
M = M+1 

BCPOS(I)=M 
POS( I)=0 
ENDIF 
C 

WRITEdl. 100) POS(I) 

WRITE(12, 100) BCPOS(I) 

C 

90 CONTINUE 
CLOSE(ll) 

CLOSE( 12) 

100 FORMAT(Il) 

ENDIF 

C 

C» 

C* 

c» INPUT CALCULATION ADJUSTMENT ITEMS. N IS THE NUMBER OF 
C* INPUT PARAMETERS IN THIS SUBMODULE. OBSERVE THAT SUBSEQUENT TO 
C* THE DEVELOPMENT OF THIS MODULE, PARAMETER NUMBER 3 WAS NO LONGER 
C* ACCESSED FROM ANY OTHER PROGRAM, AND THEREFORE SERVES ONLY AS A 
Ct DUMMY PARAMETER. THE REMAINING PARAMETERS HAVE THE FOLLOWING 
Ct MEANINGS 
Ct 

Ct 1. EPSILON REFERS TO THE DISTANCE BETWEEN THE RING LOAD 

C* AND THE BOUNDARY. NOTE THAT IN SOME INSTANCES. LETTING 
C* EPSILON EQUAL TO ZERO CAUSES A DIVIDE BY ZERO TYPE ERROR. 

Ct 

Ct 2. THE INCREMENTS FOR R REFER TO THE STATIONS ACROSS 

Ct THE PLATE. THIS APPLIES TO THE NUMBER OF INTEGRATION 
C* INCREMENTS FOR DISTRIBUTED LOADS, AS WELL AS THE NUMBER 
C* OF DATA POINTS FOR THE RESULTS. 

C* 

C* 3. THIS PARAMETER IS NO LONGER USED 

C* 

C* 4. THIS REFERS TO THE EXTENT BY WHICH THE GREENS 

Ct FUNCTION PLATE GEOMETRY OVERHANGS THE ACTUAL PLATE. 

Ct FOR EXAMPLE, A VALUE OF 0.5 MEANS THAT THE INNER RADIUS OF 
Ct THE RING LOADED PLATE DESCRIBED BY THE GREEN’S FUNCTION 
Ct OVERHANGS THE ACTUAL PLATE BY 0.5 UNITS. 

Ct 
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C* D. THE LOAD INCREMENT IS THE PERCENTAGE OF THE 

C* ORIGINALLY SPECIFIED LOAD MAGNITUDE THAT IS APPLIED DURING 
C* EACH LOAD INCREMENT AS THE PLATE IS LOADED IN THE PLASTIC 
C* RANGE 
C* 

C* 6. THIS IS THE NUMBER OF INTEGRATION INCREMENTS USED TO 

C* DETERMINE THE PLASTIC MOMENT ACROSS THE PLATE. 

C* 

C* 7. RMS IS THE ROOT MEAN SQUARE OF THE PLASTIC MOMENT 

C* INCREMENTS FOR EVERY POINT IN THE PLATE WHERE THE PLATE IS 
C* NO LONGER ELASTIC. THIS IS COMPARED WITH THE RMS VALUE FOR 
C> THE PREVIOUS INCREMENT, TO SEE WHETHER THE SOLUTION HAS 
C* SUFFICIENTLY CONVERGED. 1% CONVERGENCE MEANS THESE NUMBERS 
C* AGREE TO WITHIN LESS THAN IX 
C* 

ct 

c 

N = 7 

LADJ(l)=*EPSILON FOR RO = * 

LADJ(2)=* INCREMENTS FOR R = * 

LADJ(3)=* WIDTH INTG INCH = * 

LADJ(4)=* BOUNDS FOR A AND B = ’ 

LADJ(5)='LOAD INCREMENT = * 

LADJ(6)=*DEPTH INTG INCR = * 

LADJ(7)=*% RMS CONVERGENCE = * 

C 

WRITE(*, 102) 

102 FORMAT(//, 6X, ’WOULD YOU LIKE TO SEE THE CURRENT ADJUSTABLE’ 

+ ,/,2X, ’PARAMETER SETTINGS? ’\) 

READ(*,30) RESP 
IF (RESP .EQ. ’Y’) THEN 

OPEN ( 1 6 , FI LE= ’ AD J . D AT * , STATUS^ ’ OLD ’ ) 

WRITE(»,30) ’ ’ 

DO 108 1 = 1. N 

READ( 16, 70) ADJ( I) 

WRITE(*,106) LADJ(I) , ADJ(I) 

106 FORMAT ( 2X, A25, F20. 10) 

108 CONTINUE 

CLOSE( 16) 

ENDIF 

C 

WRITE(*,110) 

110 FORMAT(//,6X, ’DO YOU WISH TO INPUT OR CHANGE ANY COMPUTATIONAL’, 
+ /,2X, ’ADJUSTMENT ITEMS (Y/N) ? ’,\) 

READ(*,30) RESP 
C 

IF (RESP . EQ. ’ Y* ) THEN 
C 

WRITE(», 120) 

120 FORMAT(/,’ INPUT THE VALUES USING DECIMAL NOTATION ....’,/) 

DO 125 1=1, N 

WRITE(>, 126) LADJd ) 

READ(>,70) ADJ{I) 

125 CONTINUE 
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126 F0HMAT(2X, A25,\) 

C 

0PEN(13.FILE=’ADJ.DAT’ , STATUS =’ NEW ’ ) 
DO 130 1=1, N 

WHITE(13,70) ADJ(I) 

130 CONTINUE 

CLOSE( 13) 



WHITE(», 135) 

135 FORMATC///,’ YOU ABE LEAVING THE INPUT MODULE’,///) 
END 
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c 

c* 

ct PHOGHAM LOAD 

C* 

C* THIS PROGRAM ALLOWS THE USER TO INPUT THE DESIRED LOADING 

C* CONDITION FOR THE PLATE. UNIFORM DISTRIBUTED LOADS, RAMPS, AND 
C* PARABOLICALLY DISTRIBUTED LOADS ARE PERMITTED, AS WELL AS RING 
C* LOADS. TO RUN, THIS PROGRAM MUST ACCESS THE FILE LOAD.DAT. 

C* 

Ct USING THE INPUT DATA, THE PROGRAM CALCULATES THE TOTAL LOAD 

Ct PER UNIT LENGTH THAT ACTS ON A RING OF THE PLATE SURFACE OR A 
Ct WIDTH DETERMINED BY ADJUSTABLE PARAMETER #3 OF THE INPUT 
C* PROGRAM 
Ct 

Cttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 



Ct 

c 

PROGRAM LOAD 
C 

IMPLICIT REAL*8 (A-H,0-Z) 

CHARACTER*! RESP 
CHARACTER*30 LBL(4) 

REAL*8 LD(50) 

INTEGER N 
C 

LBL(1)=*MAGNITUDE AT OUTER RADIUS = * 

LBL(2)=*DISTANCE FROM PLATE CENTER = * 

' LBL(3)= ’MAGNITUDE AT INNER RADIUS = ’ 

LBL(4) =’ DISTANCE FROM PLATE CENTER = ’ 

C 

OPEN(6,FILE=’LOAD.DAT* , STATUS =* OLD * ) 

READ(6,90) LD 
N=INT(6.0+2.0*LD(6) ) 

CLOSE(6) 

C 

WRITE(*, 10) 

10 FORMAT(//,’ DO YOU WISH TO SEE THE LOADING CONDITIONS’, 
+ /,’ (Y/N) ? ’,\) 

READ(»,60) RESP 
IF (RESP .EO. ’Y’) THEN 
WHITE(», 15) 

15 FORMATC//,’ THE EXISTING LOADING CONDITIONS FOLLOW . 

DO 40 1=1,4 





WRITE(*,30) 


LBL(I) ,LD(I) 


30 


F0RMAT(2X, A30, F20. 10) 


40 


CONTINUE 
WRITER*, 60) ’ ’ 






IF (LD(5) .EQ. 


1.0) THEN 




WRITE(*,80) 

ELSE 


’LOADING IS LINEAR’ 




WRITE(*,80^ 

ENDIf 

WRITE(*,60) ’ ’ 

J = 5 


’LOADING IS PARABOLIC 



//) 
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43 

45 

C 

46 
C 

50 

60 

C 

C 

70 

75 

80 

90 

C 

100 



C 

110 

120 

130 



DO 45 I=7,6+INT(LD(6)) 

J = J + 2 

WRITE(*,60) ’ ’ 

WRITB(*,43) 'MAGNITUDE OF RING LOAD', 1-6, ' = *,LD(J) 

WRITE(*,43) 'LOCATION OF RING LOAD’, 1-6, ' = ’,LD(J+1) 

FOHMAT(2X, A24, I3,A3,F20. 10) 

CONTINUE 

ENDIF 

WRITB(»,46) 

FORMAT(//,2X, 'DO YOU WISH TO CHANGE THE LOADING CONDITIONS', 

/, ' (Y/N) ? ' ,\) 

READ(*,60) RESP 

IF (RESP ,EQ. 'Y') THEN 
WRITB(*,50) 

FORMAT(//,2X, 'DO YOU WISH TO CHANGE THE DISTRIBUTED LOAD', 

+ /.* (Y/N) ? '\) 

READ(t,60) RESP 
FORMAT( Al) 

IF (RESP .EQ. 'Y') THEN 
WRITB(*,70) 

FORMAT(//,2X, 'SPECIFY THE LOAD PROFILE USING DECIMAL', 

+ /, ' NOTATION ... ' ,//) 

DO 75 1=1,4 

WRITE(#,80) LBL(I) 

READ(»,90) LD(I) 

CONTINUE 

FORMAT(2X,A30,\) 

FORMAT(F20. 10) 

WRITE(*, 100) 

FORMATC//, 2X, ' IS THE LOAD LINEAR OR PARABOLIC (L/P) ? ',\) 
READ(*,60) RESP 
IF (RESP .BQ. 'L') THEN 
LD(5) = 1.0 
ELSE 

LD(5) = 2.0 
ENDIF 
ENDIF 

WRITE(«, 110) 

FORMAT(//,2X, 'DO YOU WISH TO MODIFY ANY RING LOADS (Y/N) ? ',\) 
READ(»,60) RESP 
IF (RESP .EQ. 'Y') THEN 
WRITE(*, 120) 

FORMAT(//,2X, 'THE NUMBER OF RING LOADS (INTEGER) = ’\) 

READ(*,130; N 

F0RMAT(I3) 

LD(6)=REAL(N) 

K = 5 

DO 200 1=1 , N 
K = K•^2 
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WRITE(»»60) * * 

WRITE(*,210) 'MAGNITUDE OF RING LOAD ' , I , »= ’ 

READ(*,90) LD(K) 

WRITE(»,210) 'LOCATION OF RING LOAD ',1, '= ' 

READ(*,90) LD(K+1) 

200 CONTINUE 

210 FORMAT(2X, A23, 13. A2,\) 

ENDIF 

C 

OPEN(7, FILE =’ LOAD. DAT' , STATUS= ' NEW ' ) 

DO 219 1=1,50 
WRITE(7,90) LD(I) 

219 CONTINUE 
ENDIF 

C 

WRITEC*, 220) 

220 FORMATC//, 6X, ' YOU ARE LEAVING THE LOADING MODULE',///) 
END 
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C 

Ct 

c» 

C* CALL THE SUBROUTINE YLD, WHICH WILL RETURN A LOAD MATRIX 

C* SUCH THAT YIELDING HAS JUST STARTED TO OCCUR IN THE PLATE 
C* 

c» 

c 

CALL YLD(KK, IMAX, LDMATl) 

C 

C* 

c* 

C* REPEAT THE PROCESS, AND THEN CALL UPDATE TO REVISE OUR STATE 

C* MATRIX TO REFLECT THE CONDITION OF THE PLATE AT THE ONSET OF 
C* YIELDING 

Ct 

Cttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

CALL AMATR(BFLAG1,BFLAG2, BC,BCPOS,POS, LDMATl, LHS, AMAT) 

CALL GAUSS( LHS, AMAT, X) 

CALL BCM(POS,BC,X,BCMAT) 

CALL SOLVER(LDMATl ,BCMAT) 

CALL UPDATE (STTl , IMAX , OUT , LDMATl , LL, LDMATO) 

CALL RESUL(M, STTl , COUNT, LDMATO) 

C 

Ct 

Cttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

ct NOW, INPUT AN INCREMENTAL LOAD, AND COMPUTE THE RESULTING 

Ct MOMENTS, DEFLECTIONS BTC, USING AS AN INPUT THE PLASTIC MOMENTS 
Ct (IF ANY) AND OTHER DATA FROM THE PREVIOUS LOAD CASE 

Ct 

Cttttttttttvtttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

c 

DO 110 1=1, NN 

LDMAT 1 ( I ) = AD J ( 5 ) t LDMATO ( I ) 

110 CONTINUE 
C 

116 CONTINUE 

WHITE (t, 130) STT(IMAX,9) ,STT(IMAX, 10) 

130 FORMAT(5X, *MRP IN = ’,F20.10,’ MTP IN = *,F20.10) 

CALL AMATR(BFLAGl,BFLAG2,BC,BCPOS,POS, LDMATl, LHS, AMAT) 

CALL GAUSS(LHS,AMAT,X) 

CALL BCM(POS,BC,X,BCMAT) 

CALL S0LVER(LDMAT1 ,BCMAT) 
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C 

C* 

c*******t***t*****t****t**t ************************** t***********t***** 
c* 

C* GO INTO THE SUBROUTINE PLAST TO CALCULATE THE INCREMENTAL PLASTIC 
C* MOMENTS AT EACH STATION ACROSS THE PLATE WHERE YIELDING HAS 
C* OCCURRED. THE PLASTIC MOMENTS THAT ARE RETURNED BY PLAST ARE 
C* SQUARED AND ADDED TO THE TOTAL IN THE VARIABLE RMS. RMS IS 
C* COMPARED TO THE RMS VALUE FROM THE PREVIOUS ITERATION (SAVED AS 
C* RMSl) 

C* 

C* ******************** ************************************************* 
C* 

C 

RMSl^RMS 
RMS=0. 0 

IZ=NINT(ADJ(6) ) 

DO 200 1=1, NN 

STTld, 1)=STT(I,9) 

STT1(I,2)=STT(I, 10) 

CALL PLAST(I) 

IF (SET(I,IZ) .GE. GEOM(6)) THEN 

RMS = RMS + STT(I, 9)**2. 0 + STT( I, 10)**2. 0 
ENDIF 
200 CONTINUE 
C 

WRITE(*,810) RMS 
810 FORMAT(5X, * RMS = *,F20.10) 

* 

C* THIS WAS THROWN IN TO AVOID A DIVIDE BY ZERO, UNTIL THE PROCESS 
C* HAS BEEN "PRIMED” 

C* 

C 

IF (RMS .EQ. 0.0) THEN 
RMS =1.0 
ENDIF 
C 

IF (ABS(SQRT(RMS)-SQRT(RMS1))/SQRT(RMS)*100.0 .GT. ADJ(7)) THEN 
DO 815 1=1, NN 

STT(I, 9) = (STT(I, 9)+STTl(I, 1 ) )/2. 0 
STT(I,10)=(STT(I, 10)+STT1( I, 2) )/2.0 
815 CONTINUE 
GOTO 116 
ENDIF 
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c 

c* 

Ct 

ct CALL THE SUBROUTINE UPDATE, WHICH WILL UPDATE OUR STATE 

C* MATRIX AND OUR LOADING MATRIX TO REFLECT THE CONDITION OF THE 
C* PLATE AT THE END OF THE LOADING INCREMENT. TEE PARAMETER 
C» COUNT KEEPS TRACK OF LOAD INCREMENTS FOR PURPOSES OF DISPLAYING 
Ct RESULTS. WHEN COUNT=10, RESULTS ARE DISPLAYED AND WRITTEN TO 
C» FILE. (IE EVERY 10 LOAD INCREMENTS) 

C* 

Ctttttxttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

c 

CALL UPDATE (STTl , IMAX , OUT , LDMAT 1 , LL, LDMATO) 

CALL RESULCM, STTl, COUNT, LDMATO) 

IF (COUNT .EQ. 10.0) THEN 
C0UNT=0, 0 
BNDIF 

COUNT=COUNT+l . 0 
M=M+1 

WRITE(*,117) M 

117 FORMAT(5X, * LOAD INCH = *,I3) 

C 

C* 

C**********»C*****C***********»*<:**»»»*»***»**»»*»*»**»»*»***C***»** 

C* 

C* CHECK TO SEE IF EXCEEDING ULTIMATE STRAIN. IF SO, KICK 

C* OUT OF THE PROCEDURE. IF NOT, WE SHOULD INCREMENT ON UP AGAIN 
C* BY USING AS OUR INPUT THE INCREMENTAL LOAD MATRIX. 

Ct 

Cttttt*tttttttttttttttttttt*tttttttttttttttttttttttttttttttttttttttt 

Ct 

IF (OUT .NE. 1) THEN 
GOTO 116 
ENDIF 
C 

CALL RESUL(M, STTl , COUNT, LDMATO) 

C 

CLOSE(15) 

C 

END 
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t 

t SUBROUTINE INTEG 

X 

* THIS SUBROUTINE IS CALLED BY THE PROGRAM MAIN TO CREATE THE 

* ARRAY LDMAT. THIS ARRAY CONTAINS THE VALUES OF THE TOTAL LOAD 

* FOR EACH SEGMENT ALONG THE PLATE RADIUS. IT IS CONSTRUCTED 

* BY USING THE AVERAGE OF THE TWO ENDPOINTS FOR DISTRIBUTED 

* LOADS AND ADDING TO THIS VALUE THE MAGNITUDE OF ANY RING LOADS 

t THAT ARE SITUATED WITHIN THE SEGMENT. 

t 

t 

SUBROUTINE INTEG ( LD . LDMAT ) 

IMPLICIT REAL»8 (A-H,0-Z) 

REAL*8 LDMAT(50) , LD(50) ,NU 
INTEGER N 

INCLUDE: ’COMMON. FOR’ 

A=GEOM(l) 

B=GEOM(2) 

NU=GEOM(3) 

D=GEOM(4) 

K=NINT(ADJ(2)) 

DEL=( A-B)/ADJ(2) 

R = B 

N=NINT(6.0-^2.0*LD(6)) 

t 

ct 

C* CHECK TO SEE IF THERE IS A DISTRIBUTED LOAD. IF NOT, PROCEED 
Ct TO THE SECTION FOR RING LOADS 
Ct 

Ctttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

c 

IF (LD(1) .NE. 0.0 .OR. LD(3) .NB. 0.0) THEN 
C 

R = B 

SLP=(L0(1)-LD(3))/(LD(2)-LD(4)) 

DO 100 1=1, K 
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c 

c* 



c* 

ct 

Ct 

ct 

c 

c 

c 



SECTION FOR LINEARLY DISTRIBUTED LOADS 

IF (LD(5) .EQ. 1.0 ) THEN 
H=R+DEL 

IF (R-LD(2) .GT. DEL/10.0) THEN 
LDMAT(I)=0.0 



ELSE IF (R .GT. LD(4)) THEN 

LDMAT(I)=(LP»(R-LD(4)-DEL/2. 0)+LD(3) )*DEL 

ELSE 

LDMAT(I)=0.0 

ENDIF 



C 
C 
C 
C 

C* 

C**»**»»**********»*»*»**********»*»******»**»**»************»******* 

c» 

ct SECTION FOR PARABOLIC LOADS (Y=X*»2) 

C* 



C* 

C 



ELSE IF (LD(5) .EQ. 2.0 ) THEN 
H=R-fDEL 

P=(LD(2)-LD(4) )*»2. 0/4. 0/(LD(l)-LD(3) ) 

IF (LD(1) .GT. LD(3)) THEN 

Y=(R-DEL/2. 0-LD(4) ) **2 . 0/4 . 0/P+LD ( 3 ) 
ELSE 

Y=-(LD(2)-H-fDEL/2. 0)*»2. 0/4. 0/P+LD(4) 
ENDIF 

IF (R-LD(2) .GT. DEL/10.0) THEN 
LDMAT(I)=0.0 

ELSE IF (R .GT. LD(4)) THEN 
LDMAT( I)=Y*DEL 

ELSE 



C 

C 



LDMAT( I)=0. 0 

ENDIF 

ENDIF 
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c 

100 CONTINUE 

ENDIF 
C 

ct 

c» 

C* SECTION FOR ADDING IN RING LOADS 

Ct 

ctttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

ct 

c 

J = 5 

IF (N .NB. 0) THEN 

DO 300 I = 7,6*^NINT(LD(6)) 

J = J + 2 

L=INT( (LD(J+1)-B)/DEL)+1 
LDMaT(L)=LDMAT( L)+LD( J) 

300 CONTINUE 

ENDIF 
C 
C 

RETURN 

END 



C 

C* 

Ctttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

ct SUBROUTINE AMATR 
C* 

C* THIS SUBROUTINE IS CALLED BY THE PROGRAM MAIN TO CREATE THE 

C» 4X4 A matrix, WHICH IS SOLVED USING GAUSS ELIMINATION TO 
C* DETERMINE THE UNKNOWN BOUNDARY CONDITIONS. 

Ct 

Cttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt*tttttt 

Ct 

c 

SUBROUTINE AMATR ( BFLAGl , BFLAG2 , BC , BCPOS , POS , LDMATl , LHS, AMAT) 

C 

IMPLICIT REAL*8 (A-H,0-2) 

CHARACTER*! BFLAG1(4), BFLAG2(4) 

REAL»8 BC(4) , LHS{4 ) . INC , NU , GRNMTl ( 4 , 6 ) , GRNMT2(4 . 6) , 

+ AMAT (4, 4) , LDMAT(4) , LDTOT(4) , LDMATl (50) ,MPRTOT(4) ,MPTTOT(4) , 

-t- STTAV(50,2) 

INTEGER POS(8) , BCPOS(8) ,M,N 
C 

$INCLUDE: ’COMMON. FOR’ 

C 

A=GEOM( 1 ) 
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B=GE0M(2) 

NU=GB0M(3) 

D=GE0M(4) 

INC=ADJ(1) 

C 

DO 101 1=1,4 
LHS(I)=0.0 

IF (I ,EQ. 1 .OR. I .EQ. 3) THEN 
RO=A-INC 
ELSE 

RO=B+INC 

ENDIF 

C 

E = A 

CALL GRNFNC(R,R0,ND,D,GRNMT1) 

R = B 

CALL GRNFNC(R,RO,NU,D,GSNMT2) 

C 

C* 

ct 

ct REDUCE THE TERMS ORIGINALLY ON THE RHS TO A 4X4 MATRIX 

Ct COEFFICIENTS (AMAT) AND A 4X1 MATRIX OF NON-COEFFICIENT 
Ct TERMS (LHS) 

C* 

Czttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

c 

IF (I .EQ. 1 .OR. I .BO. 2) THEN 
K = 0 
L=0 

DO 90 J=l,4 

IF (BFLAGl(J) .NE. *U*) THEN 
K = K+1 

LHS(I)=GRNMT1(J, 1 ) *BC ( K ) *A» ( - 1 . 0 ) * * ( J- 1 ) + LHS ( I ) 

ELSE 

L = L+1 

AMATd, L)=GRNMT1(J, 1 ) AA* (-1 . 0) ( J) 

ENDIF 

73 FORMAT(2X, F10.5) 

90 CONTINUE 

C 

DO 100 J=l,4 

IF (BFLAG2(J) .NB. *U’) THEN 
K = K+1 

LHS(I)=GRNMT2( J, 1 ) *BC (K ) *B * ( * 1 ) ** J-^LHS ( I ) 

ELSE 

L = L+1 

AMATd , L)=GRNMT2( J, 1)*B*(-1)**(J-1) 

ENDIF 

100 CONTINUE 



ELSE 
K = 0 
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L = 0 

DO 110 J=l,4 

IF (BFLAGl(J) .NE. ’U*) THEN 
K = K+1 

LHS(I)=GRNMT1(J,2)»BC(K)*A»(-1. 0)»*(J~1)+LHS(I) 
ELSE 

L = L+1 

AMAT(I,L)=GRNMT1(J,2)»A*(-1.0)*»(J) 

ENDIF 

110 CONTINUE 

C 

DO 115 J=l,4 

IF (BFLAG2(J) . NE . *U*) THEN 

K = K*^1 

LHS(I)=GRNMT2(J, 2) >BC(K)<B »(-!)»<( J)+LHS(I) 

ELSE 

L=L+1 

AMAT( I ,L)=GRNMT2( J, 2)»B»(~1) ♦♦(J-1) 

ENDIF 

115 CONTINUE 

ENDIF 
C 

C* 



c* 



ct MOVE THE TERMS ORIGINALLY ON THE LHS TO THE RHS , FORM A NEW AMAT 
C* MATRIX AND A NEW LHS MATRIX WHICH CORRESPOND TO THE A AND B MATRICES 
Ct OF THE SYSTEM AX=B 
C* 



Ct 



C 



C 



IF (BFLAGK4) . NE . ’U’ .AND. I .EQ. 1) THEN 

LHS( 1) =LHS( l)+RO»BC(BCPOS(4) ) 

ELSE IF (BFLAGK4) .EQ. *U’ .AND. I .EQ. 1) THEN 
AMAT(I,POS(4) )=AMAT(I,POS(4) )-RO 
ENDIF 



IF (BFLAG2(4) .NE. *U* .AND. I .EQ. 2) THEN 
LHS(2) =LHS(2)+RO*BC(BCPOS(8) ) 

ELSE IF (BFLAG2(4) .EQ. ’U* .AND. I .EQ. 2) THEN 
AMAT(I,POS(8) )=AMAT(I,POS(8) )~RO 
ENDIF 



IF (BFLAG1(3) .NE. *U* .AND. I .EQ. 3) THEN 

LHS ( 3 ) = LHS ( 3 ) -RO» B C ( BC POS ( 3 ) ) 

ENDIF 

IF (BFLAGK4) . NE . *U* .AND. I .EQ. 3) THEN 

LHS C 3 ) = LHS ( 3 ) + BC ( BCPOS ( 4 ) ) 

ENDIF 

IF (BFLAGK3) .EQ. ’U* .AND. I .EQ. 3) THEN 

AMAT( I , POS( 3) ) =AMAT( I , POS(3) )-^RO 
ENDIF 

IF (BFLAGK4) .EQ. ’U* .AND. I .EQ. 3) THEN 

AMAT( I , POS(4 ) ) =AMATf I , POS ( 4 ) )-l . 0 
ENDIF 
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If (BFLAG2(3) .NE. ’D’ .AND. I , 
IHS(4)=LHS(4)-R0»BC(BCP0S(7) ) 
BNDIF 


.EQ. 


4) 


THEN 


IF (BFLAG2(4) . NE . *U’ .AND. I , 

LBS ( 4 ) = LBS ( 4 ) +BC ( BC POS ( 8 ) ) 
BNDIF 


. EQ. 


4) 


THEN 


IF (BFLAG2(3) .EQ. *U’ .AND. I .EQ. 

AMAT ( I , POS ( 7 ) ) = AMAT ( I , POS ( 7 ) ) -i-RO 
ENDIF 


4) 


THEN 


IF (BFLAG2(4) .EO. ’O’ .AND. I , 


.EQ. 


4) 


THEN 



AMAT ( I , POS ( 8 ) ) = AMAT (I,POS(8))-1.0 
KNDIF 



C 

101 CONTINUE 
C 

C* 

C* 

Ct FINALLY, ADJUST THE LBS TO ACCOUNT FOR THE EXTERNAL LOAD 

C* AND PLASTICITY BY CALCULATING AND SUBTRACTING APPROPRIATE 
C* TERMS. 

C* 

C*tt**ttttttttt*tt*tttt*tttt*ttt***ttt*tttttt*t**t*t*tt**ttt*t*** 

C* 

C 

K=NINT(ADJ(2) ) 

DEL=(A-B)/ADJ(2) 

C 

DO 500 1=1, K 

R=B-DEL/2. O+REALC I)*DEL 
R0=A“ADJ(1) 

STTAVd, 1) = (STT(I, 9)+STT(I + l,9))/2.0 
STTAV( I, 2)=(STT( I, 10)+STT( I+l , 10) )/2. 0 
CALL GRNFNC(R,R0,NU,D,GRNMT1) 

LDTOT( 1)=LDT0T( 1)+LDMAT1( I) *GRNMT1( 1, 1)*R 
MPRTOT(l)=MPRTOT(l)+STTAV(I, 1)*(-GRNMT1(3, 1)/D+ 
NU/R*GRNMT1 (2, 1) )*DEL*R 

MPTTOTC 1)=MPTT0T(1)+STTAV(I,2)*(-GRNMT1(2, 1)*DEL) 
LDT0T(3)=LDT0T(3)+LDMAT1( I)»GRNMT1( 1,2)*R 
MPRTOT(3)=MPRTOT(3)+STTAV( I, l)*(-GRNMTl (3. 2)/D+ 
NU/R*GRNMT1(2,2) )*DEL»R 

MPTT0T(3)=MPTT0T(3)+STTAV(I, 2)*(“GRNMT1 (2,2) *DEL) 
RO=B+ADJ( 1) 

CALL GRNFNC(R,R0,NU,D,GRNMT1) 

LDTOT(2) =LDT0T(2)+LDMAT1 ( I)*GRNMT1 ( 1 , 1 )*R 
MPRT0T(2)=MPRT0T(2)+STTAV(I, 1 )»(-GRNMT1 (3 , 1 )/D+ 
NU/R*GRNMT1(2, 1))»DEL»R 

MPTTOT(2)=MPTTOT(2)-^STTAV( 1 , 2)<(-GRNMTl (2, 1 )»DEL) 
LDT0T(4)=LDT0T(4)+LDMAT1 ( I ) *GRNMT 1 ( 1 , 2 ) <R 
MPRTOT(4)=MPRTOT(4)+STTAV( I, 1 ) * ( -GRNMTl ( 3 , 2)/D+ 

- NU/R*GRNMT1 (2, 2) ) *DEL»R 

MPTT0T(4)=MPTT0T(4)+STTAV( I, 2)*(-GRNMTl (2, 2) *DEL) 

500 CONTINUE 
C 
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DO 105 1=1,4 

LHS ( I ) =LHS ( I ) “LDTOT( I ) +MPRTOT( I ) +MPTTOT( I ) 
105 CONTINUE 
C 

DO 1000 1=1,4 
LDTOT(I)=0.0 
MPRTOT(I)=0.0 
MPTTOT(I)=0.0 
1000 CONTINUE 
C 

RETURN 

END 



C 

C* 

c* 

ct SUBROUTINE GRNFNC 

c» 

Ct THIS SUBROUTINE CALCULATES THE VALUES OF DEFLECTION, SLOPE, 

Ct MOMENT AND SHEAR AND ASSOCIATED DERIVATIVES FOR THE GREEN’S 
Ct FUNCTION, AND RETURNS THESE VALUES TO THE MAIN PROGRAM. TO 
Ct INDICATE DERIVATIVES OF PARAMETERS, THE NUMBER 1, 2, OR 3 IS 
C* TACKED ON TO THE END OF THE VARIABLE NAME. HENCE, THE FIRST 
Ct derivative of L9 WRT RO is LSI 
c* 
ct 

Cttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

ct 

c 

SUBROUTINE GRNFNC ( R , RO , NU , D , GRNMAT ) 

C 

IMPLICIT REAL»8 (A-H,0-Z) 

REAL<8 NU, GRNMAT (4, 6) , L3 , L3 1 , L32 , L33 , LS , L9 1 , L92 , L93 
SINCLUDE: ’COMMON. FOR’ 

C 

W=1.0 

A=GEOM(l)+ADJ(4) 

B=GEOM(2)-ADJ(4) 

THB=0.0 
MB=0. 0 
QB=0. 0 
C 

IF (R .GT. RO) THEN 
RRO=1.0 
ELSE 

RRO=0. 0 
ENDIF 
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Ct*****ttt:*ttt *************** t**t**********t*t*ttt*tt***t*t***t*t*tt* 

c* 

c* CALCOLATB PLATE CONSTANTS AND PLATE FONCTIONS, AND THKIH 

C» DEHIVATIVBS 
C* 

CX»««»**t»«X«»«»»«A»X«t««**A«**«»*****X***««**********«*SA«t »«*«««»«$ 

CX 

c 

L3=RO/4. 0/AX( ( (HO/A)XX2. 0+1.0)XDLOG(A/HO)^-(HO/A)XX2.0-1. 0) 
L31=1.0/4.0/AX(DLOG(A/RO)x(3.0xhoxx2.0/AXX2.0+1.0)+2.0X 
+ (EOXX2.0/AXX2. 0-1.0)) 

L32=1.0/4.0/AX(6.0xRO/AXX2.0XDLOG(A/RO)+RO/AXx2.0-1.0/RO) 
L33=1.0/4.0/AX(6.0/AXX2.0X(DLOG(A/RO)-1.0)+1.0/AXX2.0+1. O/RO 
+ XX2.0) 

C 

L9=RO/AX((1.0+NU)/2.0XDLOG(A/HO)+(1.0-NU)/4.0X(1.0-(RO/A)XX2.0)) 

L91=1.0/4.0/AX(2.0X(1.0+NU)XDLOG(A/RO)-(3.0XNU+1.0)-3.0X(1.0-NU) 

+ XROXX2. 0/AXX2. 0) 

L92=-(1.0+ND)/2.0/A/RO-3.0XROX(1.0-NU)/2.0/AXX3. 0 
L93=(1.0+ND)/2.0/A/ROXX2.0-3.0X(1.0-NU)/2.0/AXX3.0 
C 

C1=(1.0+ND)/2.0XB/AXDLOG(A/B)+(1.0-NU)/4.0x(A/B-B/A) 

C7=1.0/2.0x(1.0-NUXx2.0)X(A/B-B/A) 

C 

G3=RO/4.0/RX(((RO/R)XX2.0+1.0)xdLOG(R/RO)i-(RO/R)Xx2.0-1.0)xRRO 

G31=1.0/4.0/RX(3.0XROXX2.0/RXX2.0xDLOG(R/RO)+2.0XROXx2.0/RXx2.0+ 

+ DLOG(R/RO)-2.0)xRRO 

G32=1.0/4.0/RX(6.0xrO/RXX2.0XDLOG(R/RO)+RO/RXX2.0-1.0/RO)XRRO 
G33 = 1.0/4.0/RX(6.0/RXX2.0X(DLOG(R/RO)-1.0)+1. 0/RXX2.0 + 1.0/ROXX2. 0) 
+ XRRO 
C 

G6=RO/4. 0/RX( (RO/R) XX2. 0-1. 0+2. OxDLOG(R/RO) ) XRRO 
G61=1.0/4.0/Rx(3.0x(ROxx2.0/RXX2.0-1.0)+2. OXDLOG(R/RO))XRRO 
G62=(3. OXRO/2. 0/RXX3. 0-1. 0/2. O/RO/R) XRRO 
G63=(3.0/2.0 /Hxx3.0+1.0/2.0/H/ROXX2.0)XRRO 
C 

G9 = RO/RX( (1.0 + NU)/2.0XDLOG(R/HO) + (1.0-NU)/4 . 0X( 1 . 0- ( RO/R) Xx2 . 0) ) 

+ XRRO 

G91=1.0/4.0/RX(2.0X(1.0+NU)X(DLOG(H/HO)-1.0)+(1.0-ND)-3.0X 
+ (1.0-NU)XROXX2.0/RXX2.0)XHRO 
G92=-1.0/4.0/RX(2.0X(1.0+NU)/RO+6.0XRO/RXX2.0X(1.0-NU))XRRO 
G93=1.0/4.0/RX(2.0X(1.0+NU)/ROXX2.0-6.0/RXX2.0X(1.0-NU))XRRO 
C 

Fl = ( 1. 0 + NU)/2. OXB/RXDLOG(R/B) + (1. 0-NU) /4. Ox(R/B-B/R) 
F2=1.0/4.0x(1.0-(B/R)xx2.0x(l. 0+2 . OXDLOG ( R/B ) ) ) 

F3 = B/4. 0/RX( ( (B/R)xx2. 0+1. 0 )xDLOG(R/B) + (B/R)xx2. 0-1 .0) 

F4=1.0/2. 0x( ( 1 . 0 +NU)xb/R+( 1 . 0-NU)XR/B) 

F5 =1.0/2. OX (1.0-(B/R)xx2.0) 

F6=B/4.0/Rx( (B/R)XX2. 0- 1 . 0+2 . OXDLOG { R/B ) ) 

F7=1.0/2.0X(1. 0-NUXx2.0)X(R/B-B/R) 

F8= 1 . 0 / 2 . OX ( 1.0+NU+(1.0-NU)x(B/R)xx2.0) 

F9 = B/RX( (1 . 0 + NU)/2. OXDLOG(R/B)-r( 1. 0-NU)/4.0x( 1 .0-(B/R)XX2.0)) 
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c 

C* 



ct 



ct CALCULATE THE DEFLECTION, SLOPE, MOMENT AND SHEAR, AND 

C* ASSOCIATED DERIVATIVES, AND STORE IN THE MATRIX GRNMAT. 

C* THIS MATRIX IS A 4 BY 4 MATRIX, WHERE ROW POSITIONS ARE AS 
C* FOLLOWS 



Ct 

Ct 

Ct 

Ct 

Ct 

Ct 

Ct 



1 DEFLECTION 

2 SLOPE 

3 MOMENT 

4 SHEAR 

AND COLUMN POSITIONS ARE AS FOLLOWS 



Ct 



ct 


1 


GfiEEN’S FUNCTION (NO DERI 


VATIVES) 






ct 


2 


1ST 


DERIVATIVE 


OF 


GREEN* 


S 


FUNCTION 


WRT 


RO 


ct 


3 


2ND 


DERIVATIVE 


OF 


GREEN* 


S 


FUNCTION 


WRT 


RO 


ct 


4 


3RD 


DERIVATIVE 


OF 


GREEN* 


S 


FUNCTION 


WRT 


RO 



Ct 

Ct FOR EXAMPLE, THE VALUE RETURNED BY GRNMAT(2,3) 

C* THE SECOND DERIVATIVE OF THE SLOPE FOR THE GREEN 
Ct WRT RO 
C* 



REPRESENTS 
S FUNCTION 



Ctttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 



Ct 

c 



YB=-W»A»*3 
THB=W*A*»2 
GRNMAT ( 1 , 1 
GRNMAT(2, 1 
GRNMAT ( 3, 1 
GRNMAT(4, 1 



.0/D»(Cl»L9/C7-L3) 

. 0/D/C7*L9 

)^-(YB+THB»R*Fl-W*R**3. 0/D*G3) 
)=THB*F4-W*R»*2. 0/D*G6 
)=THB*D/R»F7-W*R*G9 
)=W*RO/R*RRO 



YBl=-W»At»3. 

THB1=W»A**2. 

GRNMAT(1,2)= 

GRNMAT(2,2)= 

GRNMAT(3,2)= 

GRNMAT(4,2)= 



0/Dt(CKL91/C7-L31) 

0/D/C7»L91 

-(YB1+THB1*R*F1-W*R*»3. 0/D*G31) 
THB1*F4-W*R**2. 0/D*G61 
THB1*D/R»F7-W*R*G91 
W/R*RRO 



YB2=-W»A**3. 0/D*(Cl/C7»L92-L32) 

THB2=W*A*»2. 0/D/C7*L92 

GRNMAT ( 1 , 3)=-( YB2+RtFl*THB2-W*R**3. 0/D*G32) 
GRNMATC2, 3) =F4*THB2“W*R»*2. 0/D*G62 
GRNMAT(3, 3) =D*F7/R*THB2-W»R»G92 
GRNMAT(4,3)=0 



YB3=-W*A»*3. 0/D*(Cl/C7*L93-L33) 

THB3=W*A**2. 0/D/C7*L93 

GRNMAT ( 1 , 4)=-( YB3-^R*F1»THB3-W»R**3. 0/D*G3 3; 
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GBKMAT(2,4)=F4*THB3-W*B»*2. 0/D*G63 
GBNMATO, 4)=D»F7/B*THB3-W*R*G93 
GBNMAT(4.4)=0 

BBTUBN 

END 



C 

Ct 

Ct 

ct SUBBOUTINK GAUSS 
Ct 

Ct THIS SUBROUTINE USES GAUSS ELIMINATION TO DETERMINE THE 

C* UNKNOWN BOUNDARY CONDITIONS. THE MATRICES A AND B OF THE 

Ct EQUATION AX=B ARE PROVIDED BY THE SUBROUTINE AMATR. THIS 

Ct SUBROUTINE THEN RETURNS THE MATRIX X OF UNKNOWN BOUNDARY 
C* CONDITIONS. THIS SUBROUTINE WAS TAKEN FROM THE BOOK 
Ct •’NUMERICAL ANALYSIS" BY L.W. JOHNSON AND R.D. RIESS. 

Ct 

Ctttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt*t 

Ct 

c 

SUBROUTINE GAUSS(B,A,X) 

C 

IMPLICIT REAL*8 (A-H,0“Z) 

CHARACTERIZE BCLBL(8) 

REAL*8 A(4,4),B(4),X(4) ,BCMAT(8) ,BC(4) 

DIMENSION AUG(50,51) 

INTEGER P0S(8) ,BCP0S(8) ,M,N 
C 

NM1 = 3 
NP = 5 
N = 4 

SET UP THE AUGMENTED MATRIX FOR AX=B 

DO 2 1=1, N 
DO 1 J=1,N 
AUG(I, J)=A(I, J) 

1 CONTINUE 
AUG(I,NP)=B(I) 

2 CONTINUE 
C 

C THE OUTER LOOP USES ELEMENTARY ROW OPERATIONS TO TRANSFORM 
C THE AUGMENTED MATRIX TO ECHELON FORM 
r 
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*-> CO 00 <J CD 0 % 4 *^ 00 



C 



DO 8 1=1, NMl 



SEARCH FOR THE LARGEST ENTRY IN COLUMN I, ROWS I THROUGH N 

IPIVOT IS THE ROW INDEX OF THE LARGEST ENTRY 

PIVOT=0.0 
DO 3 J=I,N 
TBMP=ABS(AUG( J. I) ) 

IF(PIVOT.GE.TEMP) GO TO 3 

PIVOT=TEMP 

IPIVOT=J 

CONTINUE 

IF(PIVOT.EQ.O.O) GO TO 13 
IFdPIVOT.EQ. I) GO TO 5 

INTERCHANGE ROW I AND ROW IPIVOT 

DO 4 K=I,NP 
TEMP=AUG(I,K) 

AUG( I, K)= AUG ( IPIVOT, K) 

AUG(IPIVOT,K)=TEMP 

CONTINUE 

ZERO ENTRIES (1+ 1 ),( i+2 N , 1 ) IN THE AUGMENTED MATRIX 



IP1=I+1 

DO 7 K=IP1,N 
Q=-AUG(K,I)/AUG(I,I) 

AUG(K, I)=0.0 
DO 6 J=IP1,NP 

AUG(K, J) =Q*AUG( I, J)+AUG(K. J) 

CONTINUE 

CONTINUE 

CONTINUE 

IF(AUG(N,N).EO.O.O) GO TO 13 

BACKSOLVE TO OBTAIN A SOLUTION TO AX=B 

X(N)=AUG(N,NP)/AUG(N, N) 

DO 10 K=1,NM1 
Q=0. 0 

DO 9 J=1,K 

0=0+AUG(N-K, NP-J)»X(NP-J) 

CONTINUE 

X(N-K)=(AUG(N-K, NP)-0)/AUG(N-K,N-K) 

0 CONTINUE 

CALCULATE THE NORM OF THE RESIDUAL VECTOR, B-AX. 
SET IERR0R=1 AND RETURN 

RSO=0. 0 

DO 12 1 = 1, N 
0 = 0.0 

DO 11 J=1,N 
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Q=Q+A(I, J)»X(J) 

11 COKTINUE 
RESI=B(I)-Q 
RMAG=ABS(RBSI) 

RSQ=RSQ+RMAG*<2 

12 CONTINUE 
RNORM=SQRT(RSQ) 

IERROR=l 

ABNORMAL RETURN REDUCTION TO ECHELON FORM PRODUCES A ZERO 

ENTRY ON THE DIAGONAL. THE MATRIX A MAY BE SINGULAR. 

13 IERROR=2 
C 

RETURN 

END 



C 

Ct 

Ct 

ct SUBROUTINE BCM 
C* 

C* THIS SUBROUTINE TAKES THE KNOWN BOUNDARY CONDITIONS STORED 

Ct BY THE PROGRAM INPUT IN THE FILE BCM. DAT AS WELL AS THE 

C* ^UNKNOWN” BOUNDARY CONDITIONS RETURNED BY THE SUBROUTINE GAUSS 

Ct and CREATES THE 8X1 MATRIX OF BOUNDARY CONDITIONS BCMAT. 

C* 

Ctttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

c 

SUBROUTINE BCM ( POS , BC , X , BCMAT ) 

C 

IMPLICIT REALMS (A-H,0-Z)C 
REALC8 BC(4) , X(4) ,BCMAT(8) 

INTEGER POS(8) 

C 

K = 0 
L = 0 

DO 120 1=1,8 

IF (POS(I) .GT. 0) THEN 
K = K+1 

BCMAT(I)=X(K) 

ELSE 

L = L-^1 

BCMATC I)=BC(L) 

ENDIF 
120 CONTINUE 
RETURN 
END 
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C 

Ct 

Ctttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

ct 

ct SUBROUTINE SOLVER 
Ct 

Ct THIS SUBROUTINE SOLVES TAKES THE COMPLETE SET OF BOUNDARY 

Ct CONDITIONS SUPPLIED BY THE SUBROUTINE BCM AND SOLVES THE PLATE 
Ct BENDING PROBLEM ACROSS THE PLATE WIDTH. IMPORTANT PARAMETERS 
Ct SUCH AS DEFLECTION, SLOPE BTC ARE STORED IN THE ARRAY STT. 

Ct 

Ctttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

c 

SUBROUTINE SOLVER ( XMT , BCMaT) 

C 

IMPLICIT REAL»8 (A-H,0-Z) 

C 

REAL»8 GRNMT1C4, 6) , GRNMT2 ( 4 , 6 ) , BCMAT ( 8 ) , DT ( 50 , 4 ) , SM(50,4) ,NU, 

+ XMT(50) ,STTAV(60, 2) 

INTEGER M,N 
SINCLUDB: *C0MM0N.F0R* 

C 

Ct 

Ctttttttttttttttttttttttttttttttttttttttzttttttttttttttttttttttttttt 

Ct 

ct REZERO THE MATRICES DT AND ST 
Ct 

Cttttttttttttttttttttxtttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

c 

DO 10 1=1,50 
DO 5 J=l,4 

SM(I, J)=0.0 
DT(I, J)=0.0 
5 CONTINUE 

10 CONTINUE 
C 

A=GEOM( 1) 

B=GBOM(2) 

NU=GEOM(3) 

D=GEOM(4) 



RO = B 

DEL=(A-B)/ADJ(2) 

N=INT(ADJC2)) 
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c 

Ct 

Ct***t*tt*ttt*tttt*tttttt******t***tt**tt*ttt*t**tt****t**ttit ******* 

c* 

C* THIS LOOP CALCULATES THE LOAD TERM AND THE PLASTIC MOMENT 

C* TERM FOR EACH EQUATION. I REFERS TO THE RADIAL STATION AND J 
C* REFERS TO THE EQUATION, WHERE THE FIRST EQUATION IS FOR 
C* DEFLECTION, THE SECOND FOR DW/DR, THE THIRD FOR D^2W/DR^'2 AND 
C* THE FOURTH FOR D"3W/DR"3. THE DERIVATIVES ARE THEN USED TO 
C* CALCULATE SLOPES, MOMENTS AND SHEARS 
C* 

C* 

c 

DO 155 1=1, N+1 
C 

IF (I .EQ. 1) THEN 
RO=B+ADJ( 1) 

ENDIF 

IF (I .EQ. (N+1)) THEN 
RO=A-ADJ(l) 

ENDIF 

DO 150 M=1,NINT(ADJ(2)) 

R=B-DEL/2. 0+REAL(M)*DEL 
CALL GRNFNC(R,R0,NU,D,GRNMT1) 

STTAV(M, 1)=(STT(M, 9) +STT(M+ 1 , 9) ) /2 . 0 
STTAV(M,2) = (STT(M, 10 ) -*-STT ( M+ 1 , 10 ) ) /2 . 0 
DO 140 J=l,4 

DT( I , J)=DT( I. J)+(XMT(M)»GRNMT1( 1 , J ) *R-DE L»S TTA V ( M , 1 ) * 

+ ( -GRNMTl ( 3 , J ) /D^NU/RCGRNMTl ( 2 , J ) ) »R 

+ +DEL*STTAV(M, 2)»GRNMT1 (2, J) ) 

140 CONTINUE 

150 CONTINUE 

IF (I .EQ. 1) THEN 
RO = B-t-DEL 
ELSE 

RO=RO+DEL 
ENDIF 
155 CONTINUE 
C 
C 

C* 

C********************************************************************* 

c* 

C* AT EACH INCREMENT OF WIDTH DEL ALONG THE PLATE RADIUS, CALCULATE 
C* THE DEFLECTION, SLOPE, MOMENT AND SHEAR, WHICH ARE STORED IN 
C* STT (EACH COLUMN REPRESENTS RO , Y, TH , M, Q, Y’’, Y’’’ CORRES- 

C* PONDING TO EACH STATION I ACROSS THE RADIUS) 

C* 

C********************** ************************************ *********** 

c* 

c 

DO 200 1=1, N+1 
C 

IF (I .EQ. 1) THEN 
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R0=B+ADJ(1) 

BNDIF 

IF (I .BQ. (N+D) THEN 
R0=A-ADJ(1) 

BNDIF 

C 

R=A 

CALL GRNFNC(R,RO,NU,D,GRNMTl) 

R = B 

CALL GRNFNC(R,RO,NU,D,GRNMT2) 

DO 180 J=l,4 
DO 170 K=l,4 

SM(I, J)=A*GRNMTl(K, J)*BCMAT(K)*(-1 . 0)»*K- 
+ B»GRNMT2(K, J)^BCMAT(K + 4)»(-l . 0)»*K+SM( I, J) 

170 CONTINUE 

SM(I, J)=SM(I, J)+DT(I, J) 

180 CONTINUE 

t 

t 

t STT(I.l) 

* STT(I,2) 

* STT(I,3) 

* STT(I,4) 

* STT(I,5) 

t STT(I,6) 

* STT(I,7) 

t STT(I,8) 

t STT(I,9) 

* STT(I,10) 
t 

tttt*tttttttt*tt**tt*tttt*t****t**t*t****tttt*tttttt***tttt*t*tttttt 

t 

STT(I, l)=RO 

STT(I,2)=SM(I, 1)/R0 

STT(I . 3)=1. 0/RO»(SM(I, 2)-STT( I. 2) ) 

STT(I,7)=1.0/RO»(SM(I,3)-2.0*STT(I,3) ) 

STT(I,8)=1. 0/RO*(SM(I,4)-3.0CSTT(I,7)) 

STT( I,4)=-D»(STT(I,7)+NU/RO*STT(I,3) ) 

STT(I.5)=D*(STT(I,8)+1. 0/ROCSTT(I,7)-I.0/R0»»2.0* 

+ STT(I,3)) 

STT(I, 6)=-D*(l . 0/RO*STT(I, 3) + NU*STT(I, 7) ) 

C 

C* 

Ct THIS IS A STEP TO PREVENT THE INTERIOR POINTS FROM 

C* HAVING THE EPSILON THROWN IN. 

C* 

C 

IF (I .EQ. 1) THEN 
RO=B+DEL 
ELSE 

RO=RO+DEL 

ENDIF 



- RADIUS OF INTEREST 

- DEFLECTION 

- DW/DR 

- RADIAL MOMENT 

- SHEAR 

- TANGENTIAL MOMENT 

- DW''2/DR"2 

- DW"3/DR"3 

- RADIAL PLASTIC MOMENT 

- TANGENTIAL PLASTIC MOMENT 
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c 

200 CONTINUE 
C 

RETURN 

END 



C 

C* 

ct 

Ct SUBROUTINE YLD 
C* 

C* THIS SUBROUTINE TAKES THE PLATE BENDING SOLUTION FOR THE INPUT 
Ct LOAD INCREMENT AND CALCULATES THE LOAD AT WHICH YIELDING OF THE 
C* OUTER FIBER WILL OCCUR. THE LOAD MATRIX IS THEN INCREASED TO 
C* THIS VALUE. THIS IS ESSENTIALLY A TIME SAVER, TO AVOID STEPPING 
C* UP THE LOAD IN THE ELASTIC RANGE. 

C* 

Ct 



C 

SUBROUTINE YLD ( KK , IMAX , LDMATl ) 

C 

IMPLICIT REALMS (A-H.O-Z) 

HEAL*8 LDMATK50) ,STRMAX,STTMAX 
INTEGER KK, IMAX, IRMAX, ITMAX 
C 

$INCLUDB: ’COMMON. FOR’ 

C 

C* 

Ct FIND THE LOCATION ACROSS THE PLATE WHERE THE 
Ct ARB THE HIGHEST 
C* 

C 

STRMAX=0. 0 
STTMAX=0. 0 

DO 122 I=1,NINT(ADJ(2) )+l 

IF (ABS(STT(I,4) ) .GT. STRMAX) THEN 
STRMAX=ABS(STT(I,4)) 

IRMAX=I 

ENDIF 

IF (ABS(STT(I,6)) .GT. STTMAX) THEN 
STTMAX=ABS(STTCI,6)) 

ITMAX=I 

ENDIF 

IF (STRMAX .GT. STTMAX) THEN 
KK = 1 

IMAX=IRMAX 

ELSE 



BENDING MOMENTS 
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KK = 2 

IMAX=ITMAX 

BNDIF 

122 CONTINUE 

C 

C* 

Ct CALCULATE THE EQUIVALENT STRESS AT THE LOCATION WHERE STRESSES 
Ct ARE HIGHEST. COMPARE IT TO THE YIELD STRESS. THEN RATIO UP THE 
C* LOAD MATRIX ACCORDINGLY 
Ct 
C 

ERMAX=STT( IMAX, 7)»GEOM(5)/2. 0 
ETMAX=STT( IMAX, 3)»GEOM(5)/2. 0/STT( IMAX, 1 ) 

SRMAX=GEOM(6)/GEOM(7)/(l. 0-GEOM(3)»»2. 0)»(ERMAX+GEOM(3) »ETMAX) 
STMAX=GE0M(6)/GE0M(7)/( 1. 0-GEOM(3)»»2. 0) ♦(ETMAX+GEOM(3)»ERMAX) 
SEMAX=SQRT(SRMAX»»2. 0-SRMAX»STMAX+STMAX»»2. 0) 

R=GEOM(6)/SEMAX 

DO 126 J=1,NINT(ADJ(2)) 

LDMAT1( J) =R»LDMAT1 ( J) 

126 CONTINUE 

WRITE(»,130) LDMATl(l) 

130 FORMAT(2X, * LDMATl = *,F20.10) 

C 

RETURN 

END 



C 

Ct 

Ctttttttttttttttttttztttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

Ct SUBROUTINE PLAST 
Ct 

Ct THIS SUBROUTINE CALCULATES PLASTIC MOMENTS USING A TRAPAZOIDAL 

Ct RULE INTEGRATION SCHEME. THE SUBROUTINE WILL ONLY HANDLE 
Ct LINEAR ELASTIC-PERFECTLY PLASTIC MATERIALS. 

Ct 

Ctttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

C 

SUBROUTINE PLAST(I) 

C 

IMPLICIT REALMS (A-H,0-Z)C 
REAL»8 MRP,MTP,NU 
INTEGER IZ 

$ INCLUDE: ’COMMON. FOR* 
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c 

c* 

********** 



TOTAL RADIAL STRAIN INCREMENT 
TOTAL TANGENTIAL STRAIN INCREMENT 
ELASTIC RADIAL STRAIN INCREMENT 
ELASTIC TANGENTIAL STRAIN INCREMENT 
DEPTH INCREMENT 

RADIAL POSITION STATION NUMBER 
HALF-THICKNESS STATION NUMBER 
RADIAL STRESS INCREMENT 
TANGENTIAL STRESS INCREMENT 

MATRIX HOLDING THE TEMPORARY STRESSES AT THE 

GIVEN LOADING INCREMENT ( 1=RADI AL; 2=TANGENTIAL) 
TEMPORARY EQUIVALENT STRESS - USED FOR YIELD 
CRITERIA 

EQUIVALENT STRESS FROM PREVIOUS INCREMENT 
ARRAY CONTAINING THE TEMPORARY STRAIN INCREMENTS 
AT THE GIVEN LOADING INCREMENT (1=RADIAL;2= 
TANGENTIAL 

CONSTANT REFERRING TO THE FRACTION BEYOND YIELDING 
THAT OCCURRED DURING THE COURSE OF THE LAST LOAD 
INCREMENT 
DEPTH 



c* 




c* 


varia; 


c* 




c* 


DER 


c* 


DBT 


c* 


DERE 


c* 


DBTE 


c* 


DEL 


c* 


I 


c* 


IZ 


c* 


DSl 


c* 


DS2 


c* 


STRT 


c* 




c* 


SET 


c* 




c* 


SE 


c* 


DEPT 


c* 




c* 




c* 


R 


c* 




c* 




c* 


Z 


c* 




c****** 


***** 


c* 





STT(I,9)=0.0 
STT(I, 10)=0.0 
SSUMR=0. 0 
SSUMT=0. 0 

DEL = GEOM(5)/2. 0/(ADJ(6)-l . 0) 
2 = 0 

E=GEOM(6)/GEOM(7) 

NU=GEOM(3) 



DO 500 IZ=NINT(ADJ(6)),1,-1 
C 

C* 

c*** ************************************************ ****** *********** 
c* 

THIS SECTION CALCULATES THE STRAIN AND STRESS INCREMENTS 
CORRESPONDING TO THE GIVEN LOAD INCREMENT. IT THEN ADDS THE 
STRESS INCREMENT TO THE PREVIOUS TOTAL STRESS TO CALCULATE THE 
EQUIVALENT STRESS. THE EQUIVALENT STRESS IS COMPARED TO THE 
UNIAXIAL YIELD STRESS TO SEE IF YIELDING HAS OCCURRED. IF SO, 

THE PLASTIC STRAIN INCREMENT IS SET EQUAL TO THE TOTAL STRAIN 
INCREMENT 



C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

c**tt***ttt**t ************************************************* ****** 
c* 
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c 

DBR=-STT(I, 7)»( IZ-1)*DEL 
DBT=-STT( I ,3)»( IZ-1) *DEL/ST( I , 1) 

DERE=DER-DEPT(I,IZ,1) 

DETE=DET-DEPT(I, IZ,2) 

DSl=E/( 1. 0-NU»»2. 0)»(DERE+NU*DETE) 

DS2=E/( 1 . 0-NU»*2. 0) *(DETE+NU»DERE) 

STRT(I, IZ, 1) = STR(I, IZ, D+DSl 
STRTd, IZ,2)=STR(I, IZ.2)+DS2 

SET(I, IZ)=SQRT(STRT(I,IZ, 1 ) **2 . 0-STRT ( I , IZ , 1 ) »S TRT ( I , IZ,2) 

+ ^-STHTd, IZ,2)»*2) 

C 

C* 

Ct SEE IF YIELD CRITERIA IS EXCEEDED. IF SO, REZERO THE PLASTIC 

Ct STRAIN INCREMENTS AND THE STRESS INCREMENTS. ALSO, RETURN TO THE 
Ct MAIN PROGRAM IF THE OUTER FIBER HAS NOT BEGUN TO YIELD 
Ct 
C 

IF (SET(I,IZ) .LT. GEOM(6)) THEN 
DEPTd, IZ, 1)=0. 0 
DEPT(I,IZ,2)=0.0 
STT(I,9)=0.0 
STTd, I0)=0.0 

IF (IZ .EQ. NINT(ADJ(6) ) ) THEN 
RETURN 
ENDIF 
ELSE 

SBTd, IZ)=GEOM(6) 

DEPTd, IZ,1)=DER 
DEPTd, IZ,2)=DET 
IF (SE(I,IZ) .LT. GEOM(6)) THEN 

R=(SET(I,IZ)-GEOM(6) )/(SET(I, IZ)-SE(I, IZ) ) 

DEPTd, IZ, 1)=R*DEPT(I, IZ, D 
DEPTd, IZ,2)=R»DEPT(I, IZ,2) 

STRTd, IZ, 1)=STR(I, IZ, 1)*^DS1»(1-R) 

STRTd, IZ,2)=STR(I , IZ , 2 ) +DS2* (1-R ) 

ENDIF 

ENDIF 

C 

Z=Z+DEL 

C 

500 CONTINUE 
C 

Ct 

Cttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

Ct 

ct PERFORM A TRAPAZOIDAL RULE INTEGRATION OF THE STRAIN OVER THE 

Ct PLATE HALF THICKNESS, AND MULTIPLY BY 2 TO GIVE THE TOTAL RADIAL 
Ct AND TANGENTIAL PLASTIC MOMENT INCREMENTS (STT(I.9) AND STT(I,10)) 
C» 

Ct 



Z=-DEL/2.0 

DO 540 IZ=1,NINT(ADJ(6)“1) 
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ZrZ^DEL 

DMPfi=( (DEPTd, IZ, 1 )-^DEPT(I, IZ + 1 , 1 ) ) '•'GEOM ( 3 ) » ( DEPT (I . IZ, 2) 
'*'DEPT(I, rZ + l,2)))*Z/2.0 

DMPT=((DEPT(I, IZ,2)^DEPT(I, IZ-^1, 2))'^GKOM(3)*(DEPT(I. IZ, D 
-^DEPTd, IZ + 1. l)))*Z/2.0 
SSUMR = SSUMR-^DMPR 
SSUMT=SSUMT+DMPT 
540 CONTINUE 
C 

STTd, 9)-2.0*GEOM(6)/GEOM(7)/(l-GEOM(3)*»2. 0)»DEL*SSUMR 
STT(I , 10)=2. 0»GEOM(6)/GEOM(7)/(1-GEOM(3) ♦<2. 0)*DBL*SSUMT 

WRITE(*,600) I 

600 FORMAT(2X, dAT POS = * , 12) 

WRITE(»,700) STT(I,9),STTd,10) 

700 F0RMAT(2X,’MRP = *,F20,12,* MTP = *,F20.12) 

C 

RETURN 

END 



C 

C» 

c* 

ct SUBROUTINE UPDATE 
C* 

C* THIS SUBROUTINE TAKES THE INCREMENTAL STRESSES, STRAINS, 

C* DEFLECTIONS, MOMENTS, ETC AND ADDS THEM TO THE PREVIOUS TOTALS 
C* ONCE CONVERGENCE OF PLASTIC MOMENTS HAS OCCURRED. 

C* 

C* 

C 

SUBROUTINE UPDATE (STTl , IMAX,OUT, LDMATl , LL, LDMATO) 

C 

IMPLICIT REALMS (A-H,0-Z) 

REAL* 8 STT1(50,2) , LDMATO (50) , LDMATl (50) ,NU,DEL, E 
INTEGER NN,OUT,LL 
SINCLUDE: ’COMMON. FOR’ 

C 

NN=NINT( ADJ(2) )+l 
E=GEOM(6)/GEOM(7) 

DEL = GEOM(5)/2. 0/(ADJ(6)-l. 0) 

NU=GEOM(3) 

C 

DO 1000 1=1, NN 
DO 905 J=LL, 10 

STd , J)=STd, J) + STT(I, J) 

905 CONTINUE 
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LDMATO ( I ) =LDMAT1 ( I ) +LDMATO ( I ) 

C 

DO 980 I2=1,NINT(ADJ(6)) 

C 

C* 

C> IF THE YIELD CRITERIA HAS BEEN EXCEEDED, THEN ADD IN THE 
C» PLASTIC STRAIN INCREMENTS AND UPDATE THE TOTAL STRESS TO 
C» EQUAL THE TEMPORARY STRESS CALCULATED IN THE SUBROUTINE 
C» PLAST 
Ct 
C 

IF (SE(I.IZ) .GT. GEOM(6)) THEN 

EP(I, IZ. 1)=EP(I , IZ, 1)+DEPT(I. IZ, 1) 

EP( I . IZ , 2)=EP( I . IZ, 2)-rDEPT( I, IZ .2) 
STR(I,IZ.l)=STRT(I,IZ.l) 

STR(I, IZ,2)=STRT(I, IZ. 2) 

2 = (GEOM(5)/2.0/(ADJ(6)“1.0) )*(IZ-‘l) 

ELSE 

C 

C> 

Ct IF THE YIELD CRITERIA HAS NOT BEEN EXCEEDED, THEN UPDATE THE 
C> STRESSES AND EQUIVALENT STRESSES USING THE ELASTIC EQUATIONS 

C» 

C 

STR(I, IZ. 1)=-E/(1.0-NU<>2.0)*(ST(I.7)»(IZ-1)<DEL- 
+ EP( I. IZ, 1)+NU»(ST(I,3)*(IZ-1)»DEL/ST(I. 1) ) ) 

STR(I, IZ,2)=-E/(1.0~NU*»2.0)»(ST(I,3)*(IZ-1)*DEL/ST(I. 1) 
EP(I. IZ,2)+NU»(ST( I,7)*(IZ-1)*DEL) ) 

SE( I , IZ)=SQRT(STR( I , IZ , 1)»>2. 0-STR(I , IZ , 1)» 

+ STR(I, IZ,2)+STR(I, IZ,2)»<2) 

ENDIF 

IF (ABS(EP(I, IZ. 1) ) .GT. GEOM(9) .OR. ABS ( EP ( I , IZ . 2 ) ) 

+ .GT. GEOM(9)) THEN 

OUT=l 
ENDIF 

980 CONTINUE 

1000 CONTINUE 
C 

LL=2 

RETURN 

END 
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c 

c» 

c******************************************************************** 

c« 

C* SUBROUTINE BESUL 
C* 

C» STORE THE RESULTS AT THE END OF EVERY 10 LOAD INCREMENTS IN 

C» THE FILE OUT. DAT. ALSO SCREEN DUMP THESE RESULTS. 

C* 

Ct**tt**ttt********t***********tt**t**t******tt*tt*ttttttt**tttttt*tt 

c» 

c 

SUBROUTINE RESUL(MM, STTl .COUNT, LDMATO) 

C 

IMPLICIT REALMS (A-H.O-Z) 

REALMS STT1(50,2) , LDMAT0(50) 

INTEGER NN 

SINCLUDE: ’COMMON. FOB’ 

C 

NN=NINT(ADJ(2) )+l 
DEL=(GEOM(l)-GBOM(2))/ADJ(2) 

C 

IF (COUNT .EQ. 1.0) THEN 
WRITE(15, 100) 

WRITE(*, 100) 

100 FORMAT (/// , ’ tt*tt***t*tt*t*ttt*t**tt*****tt****t*tt**tt***tt 

+ •«*»*«»»»»*»»*»»»»»»»»»*»*»»»♦**»• ,//) 

WRITE(15, 120) MM, LDMAT0(1)/DEL 
WBITE(»,120) MM, LDMAT0(1)/DEL 
120 FORMAT(2X, ’ LOAD INCR = ’,13,’ LOAD = ’.F10.4,/) 

C 

WRITE(15, 130) 

WRITE(», 130) 

130 FORMAT(2X, ’RADPOS’ ,2X, ’DEFLECT ’,2X, ’SLOPE ’,2X,’RAD MOM ’ 

+ ,2X,’ MRP ’,2X,’ TAN MOM ’,2X,’ MTP ’,’ SHEAR ’,/) 

DO 160 1=1, NN 

STM1=ST(I,4)+ST(I,9) 

STM2=ST(I,6)+ST(I,10) 

HRITE(15, 140) ST(1, 1) , -ST( 1 , 2) , ST ( 1 , 3 ) , STMl , ST( 1 , 9 ) , 

+ STM2,ST( I, 10) ,ST(I,5) 

WBITE(»,140) ST(I,1) ,-ST(I,2) ,ST(I,3) ,STM1,ST(I,9) , 

+ STM2,ST(I,10),ST(I,5) 

140 FORMAT(2X,F5.2,2X,F7.5,2X,F7.4,2X,F10.0,2X,F5.0,2X,F10.0, 

+ 2X,F5.0,2X,F10.0) 

160 CONTINUE 

C 

WRITE(15, 163) 

WRITE(*, 163) 

163 FORMAT(/,2X, ’ _ 

+ ’ IIIIIIIII ’) 

DO 1000 1=1, NN 

IF (ST(I,9) .NE. 0.0 .OR. ST(I,10) .NB. 0.0) THEN 
WRITE(15, 170) ST(I,1) 

WRITB(*,170) ST(I,1) 

170 FORMAT(/,2X, ’RADIAL POSITION = ’ , F5 . 2 ) 
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WHITB(15,165) 

WRITB(*, 165) 

165 rOBMAT(/,2X, ' DEPTH ’.2X,* EBP ’.2X,’ SBAD 

+ 2X,’ BTP ’,2X,’ STAN ’,2X,* SB ',/) 

DO 200 IZ=1,NINT(ADJ(6) ) 

IF (SB(I,IZ) .GT. GB0M(6)) THEN 

Z=(GBOM(5)/2.0/(ADJ(6)-1.0))*(IZ-1) 

WRITE (15, 180) Z,EP(I,IZ, 1) ,STR(I, IZ, 1) ,BP(I, IZ,2) , 
+ STR(I,IZ,2),SE(I,IZ) 

WRITB(*, 180) Z,EP(I, IZ,1) ,STR(I, IZ, 1) ,BP(I, IZ.2) , 

+ STB(I,IZ,2),SE(I,IZ) 

ISO F0BMAT(2X,F6.4,3X, F10.8, 3X, F8. 1, 3X, no. 8,3X,F8. 1 , 

+ 3X,F8.1) 

ENDIF 

200 CONTINUE 

ENDIF 

1000 CONTINUE 
ENDIF 
RETURN 
END 



COMMON GEOM(9) , ADJ(7) ,ST(50, 10) 

COMMON /PLSCM/ EP(50,50,2) , STT(50 , 10 ) , DEPT (50 , 50 , 2) , STB (50, 50, 2), 
+ SE(50,50) ,STRT(50,50,2),SET(50,50) 

REAL*8 ADJ, ST, EP,STT,DEPT,STR,SE,STRT,SET 
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